MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
strnode.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/strnode.cc,v 1.166 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 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33 
34 #include <cerrno>
35 #include <cstring>
36 #include <sstream>
37 
38 #include "mynewmem.h"
39 #include "strnode.h"
40 #include "body.h"
41 #include "autostr.h"
42 #include "dataman.h"
43 
44 #include "matvecexp.h"
45 #include "Rot.hh"
46 
47 static const char xyz[] = "xyz";
48 
49 static const char *sdn_dof[] = {
50  "position P",
51  "momentum B"
52 };
53 static const char *sdn_eq[] = {
54  "momentum definition B",
55  "force equilibrium F"
56 };
57 static const char *sdn_initial_dof[] = {
58  "position P",
59  "velocity v"
60 };
61 static const char *sdn_initial_eq[] = {
62  "position constraint P",
63  "position constraint derivative v"
64 };
65 
66 static const char *sn_dof[] = {
67  sdn_dof[0], // "position P",
68  "incremental rotation parameter g",
69  sdn_dof[1], // "momentum B",
70  "momenta moment G"
71 };
72 static const char *sn_eq[] = {
73  sdn_eq[0], // "momentum definition B",
74  "momenta moment definition G",
75  sdn_eq[1], // "force equilibrium F",
76  "moment equilibrium M"
77 };
78 static const char *sn_modal_eq[] = {
79  "linear velocity definition v",
80  "angular velocity definition w",
81  "force equilibrium F",
82  "moment equilibrium M"
83 };
84 static const char *sn_initial_dof[] = {
85  sdn_initial_dof[0], // "position P",
86  "incremental rotation parameter g",
87  sdn_initial_dof[1], // "velocity v",
88  "angular velocity w"
89 };
90 static const char *sn_initial_eq[] = {
91  sdn_initial_eq[0], // "position constraint P",
92  "orientation constraint g",
93  sdn_initial_eq[1], // "position constraint derivative v",
94  "orientation constraint derivative w"
95 };
96 
97 
98 
99 /* StructDispNode - begin */
100 
101 /* Costruttore definitivo */
102 StructDispNode::StructDispNode(unsigned int uL,
103  const DofOwner* pDO,
104  const Vec3& X0,
105  const Vec3& V0,
106  const StructNode *pRN,
107  const RigidBodyKinematics *pRBK,
108  doublereal dPosStiff,
109  doublereal dVelStiff,
110  OrientationDescription od,
111  flag fOut)
112 : Node(uL, pDO, fOut),
113 XPrev(X0),
114 XCurr(X0),
115 VPrev(V0),
116 VCurr(V0),
117 XPPCurr(Zero3),
118 XPPPrev(Zero3),
119 pRefNode(pRN),
120 #ifdef USE_NETCDF
121 Var_X(0),
122 Var_Phi(0),
123 Var_XP(0),
124 Var_Omega(0),
125 Var_XPP(0),
126 Var_OmegaP(0),
127 #endif /* USE_NETCDF */
128 od(od),
129 dPositionStiffness(dPosStiff),
130 dVelocityStiffness(dVelStiff),
131 pRefRBK(pRBK),
132 bOutputAccels(false)
133 {
134  NO_OP;
135 }
136 
137 /* Distruttore (per ora e' banale) */
138 StructDispNode::~StructDispNode(void)
139 {
140  NO_OP;
141 }
142 
143 /* Tipo di nodo */
144 Node::Type
145 StructDispNode::GetNodeType(void) const
146 {
147  return Node::STRUCTURAL;
148 }
149 
150 /* rigid-body kinematics */
151 const RigidBodyKinematics *
152 StructDispNode::pGetRBK(void) const
153 {
154  return pRefRBK;
155 }
156 
157 const Vec3&
158 StructDispNode::GetX(void) const
159 {
160  return GetXCurr();
161 }
162 
163 const Mat3x3&
164 StructDispNode::GetR(void) const
165 {
166  return ::Eye3;
167 }
168 
169 const Vec3&
170 StructDispNode::GetV(void) const
171 {
172  return GetVCurr();
173 }
174 
175 const Vec3&
176 StructDispNode::GetW(void) const
177 {
178  return ::Zero3;
179 }
180 
181 const Vec3&
182 StructDispNode::GetXPP(void) const
183 {
184  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
185 }
186 
187 const Vec3&
188 StructDispNode::GetWP(void) const
189 {
190  return ::Zero3;
191 }
192 
193 bool
194 StructDispNode::ComputeAccelerations(bool b)
195 {
196  return false;
197 }
198 
199 std::ostream&
200 StructDispNode::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
201 {
202  integer iIndex = iGetFirstIndex();
203 
204  out
205  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
206  "position [px,py,pz]" << std::endl;
207 
208  if (bInitial) {
209  iIndex += 3;
210  out
211  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
212  "linear velocity [vx,vy,vz]" << std::endl;
213  }
214 
215  return out;
216 }
217 
218 void
219 StructDispNode::DescribeDof(std::vector<std::string>& desc, bool bInitial, int i) const
220 {
221  if (i == -1) {
222  if (bInitial) {
223  desc.resize(6);
224 
225  } else {
226  desc.resize(3);
227  }
228 
229  } else {
230  desc.resize(1);
231  }
232 
233  std::ostringstream os;
234  os << "StructDispNode(" << GetLabel() << ")";
235 
236  // always uses initial_dof[] becuase dof[]
237  // and initial_dof[] are the same up to 6
238  int iend = bInitial ? 6 : 3;
239  if (i == -1) {
240  std::string name = os.str();
241 
242  for (i = 0; i < iend; i++) {
243  os.str(name);
244  os.seekp(0, std::ios_base::end);
245  os << ": " << sdn_initial_dof[i/3] << xyz[i%3];
246  desc[i] = os.str();
247  }
248 
249  } else {
250  if (i < 0 || i >= iend) {
251  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
252  }
253  os << ": " << sdn_initial_dof[i/3] << xyz[i%3];
254  desc[0] = os.str();
255  }
256 }
257 
258 std::ostream&
259 StructDispNode::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
260 {
261  integer iIndex = iGetFirstIndex();
262 
263  if (bInitial) {
264  out
265  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
266  "position [Px,Py,Pz]" << std::endl
267  << prefix << iIndex + 7 << "->" << iIndex + 9 << ": "
268  "linear velocity [vx,vy,vz]" << std::endl;
269  } else {
270  if (dynamic_cast<const DynamicStructDispNode*>(this) != 0) {
271  iIndex += 3;
272  }
273 
274  out
275  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
276  "force equilibrium [Fx,Fy,Fz]" << std::endl;
277  }
278 
279  return out;
280 }
281 
282 void
283 StructDispNode::DescribeEq(std::vector<std::string>& desc, bool bInitial, int i) const
284 {
285  if (i == -1) {
286  if (bInitial) {
287  desc.resize(6);
288 
289  } else {
290  desc.resize(3);
291  }
292 
293  } else {
294  desc.resize(1);
295  }
296 
297  std::ostringstream os;
298  os << "StructDispNode(" << GetLabel() << ")";
299 
300  if (i == -1) {
301  std::string name(os.str());
302 
303  if (bInitial) {
304  for (i = 0; i < 6; i++) {
305  os.str(name);
306  os.seekp(0, std::ios_base::end);
307  os << ": " << sdn_initial_eq[i/3] << xyz[i%3];
308  desc[i] = os.str();
309  }
310 
311  } else {
312  for (i = 0; i < 3; i++) {
313  os.str(name);
314  os.seekp(0, std::ios_base::end);
315  os << ": " << sdn_eq[1 + i/3] << xyz[i%3];
316  desc[i] = os.str();
317  }
318  }
319 
320  } else {
321  if (bInitial) {
322  if (i < 0 || i >= 6) {
323  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
324  }
325 
326  os << ": " << sdn_initial_eq[i/3] << xyz[i%3];
327 
328  } else {
329  if (i < 0 || i >= 3) {
330  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
331  }
332 
333  os << ": " << sdn_eq[1 + i/3] << xyz[i%3];
334  }
335  desc[0] = os.str();
336  }
337 }
338 
339 /* Contributo del nodo strutturale al file di restart */
340 std::ostream&
341 StructDispNode::Restart(std::ostream& out) const
342 {
343  out << " structural displacement: " << GetLabel() << ", ";
344  if (GetStructDispNodeType() == StructDispNode::DYNAMIC) {
345  out << "dynamic";
346  } else if (GetStructDispNodeType() == StructDispNode::STATIC) {
347  out << "static";
348  }
349  out << ", reference, global, ";
350  XCurr.Write(out, ", ")
351  << ", reference, global, ",
352  VCurr.Write(out, ", ")
353  << ", assembly, "
354  << dPositionStiffness << ", "
355  << dVelocityStiffness
356  << ", scale, " << pGetDofOwner()->dGetScale() << ';' << std::endl;
357 
358  return out;
359 }
360 
361 DofOrder::Order
362 StructDispNode::GetDofType(unsigned int i) const
363 {
364  ASSERT(i >= 0 && i < iGetNumDof());
365  return DofOrder::DIFFERENTIAL;
366 }
367 
368 
369 /* Restituisce il valore del dof iDof;
370  * se differenziale, iOrder puo' essere = 1 per la derivata */
371 const doublereal&
372 StructDispNode::dGetDofValue(int iDof, int iOrder) const
373 {
374  ASSERT(iDof >= 1 && iDof <= 3);
375  ASSERT(iOrder == 0 || iOrder == 1);
376  if (iDof >= 1 && iDof <= 3) {
377  if (iOrder == 0) {
378  return XCurr(iDof);
379  } else if (iOrder == 1) {
380  return VCurr(iDof);
381  }
382  } else {
383  silent_cerr("StructDispNode(" << GetLabel() << "): "
384  "required dof " << iDof << " (order " << iOrder << ") "
385  "is not available." << std::endl);
386  throw StructDispNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
387  }
388 
389  /* dummy return value to workaround compiler complains */
390  static doublereal dmy = 0.;
391  return dmy;
392 }
393 
394 /* Restituisce il valore del dof iDof al passo precedente;
395  * se differenziale, iOrder puo' essere = 1 per la derivata */
396 const doublereal&
397 StructDispNode::dGetDofValuePrev(int iDof, int iOrder) const
398 {
399  ASSERT(iDof >= 1 && iDof <= 3);
400  ASSERT(iOrder == 0 || iOrder == 1);
401  if (iDof >= 1 && iDof <= 3) {
402  if (iOrder == 0) {
403  return XPrev(iDof);
404  } else if (iOrder == 1) {
405  return VPrev(iDof);
406  }
407  } else {
408  silent_cerr("StructDispNode(" << GetLabel() << "): "
409  "required dof " << iDof << " (order " << iOrder << ") "
410  "is not available." << std::endl);
411  throw StructDispNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
412  }
413 
414  /* dummy return value to workaround compiler complains */
415  static doublereal dmy = 0.;
416  return dmy;
417 }
418 
419 /* Setta il valore del dof iDof a dValue;
420  * se differenziale, iOrder puo' essere = 1 per la derivata */
421 void
422 StructDispNode::SetDofValue(const doublereal& dValue,
423  unsigned int iDof,
424  unsigned int iOrder /* = 0 */ )
425 {
426  ASSERT(iDof >= 1 && iDof <= 6);
427  ASSERT(iOrder == 0 || iOrder == 1);
428  if (iDof >= 1 && iDof <= 3) {
429  if (iOrder == 0) {
430  XCurr(iDof) = dValue;
431 
432  } else if (iOrder == 1) {
433  VCurr(iDof) = dValue;
434  }
435 
436  } else {
437  silent_cerr("StructDispNode(" << GetLabel() << "): "
438  "required dof " << iDof << " (order " << iOrder << ") "
439  "is not available." << std::endl);
440  throw StructDispNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
441  }
442 }
443 
444 void
445 StructDispNode::OutputPrepare(OutputHandler &OH)
446 {
447  if (bToBeOutput()) {
448 #ifdef USE_NETCDF
449  if (OH.UseNetCDF(OutputHandler::STRNODES)) {
450  ASSERT(OH.IsOpen(OutputHandler::NETCDF));
451 
452  // node
453  const char *type;
454  switch (GetStructDispNodeType()) {
455  case STATIC:
456  type = "static";
457  break;
458 
459  case DYNAMIC:
460  type = "dynamic";
461  break;
462 
463  default:
464  pedantic_cerr("StructDispNode::OutputPrepare(" << GetLabel() << "): "
465  "warning, unknown node type?" << std::endl);
466  type = "unknown";
467  break;
468  }
469 
470  std::ostringstream os;
471  os << "node.struct." << GetLabel();
472  (void)OH.CreateVar(os.str(), type);
473 
474  // node sub-data
475  os << '.';
476  std::string name = os.str();
477 
478  Var_X = OH.CreateVar<Vec3>(name + "X", "m",
479  "global position vector (X, Y, Z)");
480 
481  Var_Phi = OH.CreateRotationVar(name, "", od, "global");
482 
483  Var_XP = OH.CreateVar<Vec3>(name + "XP", "m/s",
484  "global velocity vector (v_X, v_Y, v_Z)");
485 
486  Var_Omega = OH.CreateVar<Vec3>(name + "Omega", "radian/s",
487  "global angular velocity vector (omega_X, omega_Y, omega_Z)");
488 
489  // accelerations
490  if (bOutputAccels) {
491  Var_XPP = OH.CreateVar<Vec3>(name + "XPP", "m/s^2",
492  "global acceleration vector (a_X, a_Y, a_Z)");
493 
494  Var_OmegaP = OH.CreateVar<Vec3>(name + "OmegaP", "radian/s^2",
495  "global angular acceleration vector (omegaP_X, omegaP_Y, omegaP_Z)");
496  }
497  }
498 #endif // USE_NETCDF
499  }
500 }
501 
502 /* Output del nodo strutturale (da mettere a punto) */
503 void
504 StructDispNode::Output(OutputHandler& OH) const
505 {
506  if (bToBeOutput()) {
507 #ifdef USE_NETCDF
508  if (OH.UseNetCDF(OutputHandler::STRNODES)) {
509  Var_X->put_rec(XCurr.pGetVec(), OH.GetCurrentStep());
510  switch (od) {
511  case EULER_123:
512  case EULER_313:
513  case EULER_321:
514  case ORIENTATION_VECTOR:
515  Var_Phi->put_rec(::Zero3.pGetVec(), OH.GetCurrentStep());
516  break;
517 
518  case ORIENTATION_MATRIX:
519  Var_Phi->put_rec(::Eye3.pGetMat(), OH.GetCurrentStep());
520  break;
521 
522  default:
523  /* impossible */
524  break;
525  }
526  Var_XP->put_rec(VCurr.pGetVec(), OH.GetCurrentStep());
527  Var_Omega->put_rec(::Zero3.pGetVec(), OH.GetCurrentStep());
528 
529  if (bOutputAccels) {
530  Var_XPP->put_rec(XPPCurr.pGetVec(), OH.GetCurrentStep());
531  Var_OmegaP->put_rec(::Zero3.pGetVec(), OH.GetCurrentStep());
532  }
533  }
534 #endif /* USE_NETCDF */
535 
536  if (OH.UseText(OutputHandler::STRNODES)) {
537  std::ostream& out = OH.StrNodes();
538  out
539  << std::setw(8) << GetLabel()
540  << " " << XCurr << " ";
541  switch (od) {
542  case EULER_123:
543  case EULER_313:
544  case EULER_321:
545  case ORIENTATION_VECTOR:
546  OH.StrNodes() << ::Zero3;
547  break;
548 
549  case ORIENTATION_MATRIX:
550  OH.StrNodes() << ::Eye3;
551  break;
552 
553  default:
554  /* impossible */
555  break;
556  }
557  OH.StrNodes() << " " << VCurr << " " << ::Zero3;
558 
559  if (bOutputAccels) {
560  out
561  << " " << XPPCurr
562  << " " << ::Zero3;
563  }
564  out << std::endl;
565  }
566  }
567 }
568 
569 /* Aggiorna dati in base alla soluzione */
570 void
571 StructDispNode::Update(const VectorHandler& X, const VectorHandler& XP)
572 {
573  integer iFirstIndex = iGetFirstIndex();
574 
575  XCurr = Vec3(X, iFirstIndex + 1);
576  VCurr = Vec3(XP, iFirstIndex + 1);
577 }
578 
579 
580 /* Aggiorna dati in base alla soluzione */
581 void
582 StructDispNode::DerivativesUpdate(const VectorHandler& X, const VectorHandler& XP)
583 {
584  integer iFirstIndex = iGetFirstIndex();
585 
586  /* Forza configurazione e velocita' al valore iniziale */
587  const_cast<VectorHandler &>(X).Put(iFirstIndex + 1, XCurr);
588  const_cast<VectorHandler &>(XP).Put(iFirstIndex + 1, VCurr);
589 }
590 
591 
592 /* Aggiorna dati in base alla soluzione durante l'assemblaggio iniziale */
593 void
594 StructDispNode::InitialUpdate(const VectorHandler& X)
595 {
596  integer iFirstIndex = iGetFirstIndex();
597 
598  XCurr = Vec3(X, iFirstIndex + 1);
599  VCurr = Vec3(X, iFirstIndex + 4);
600 }
601 
602 /* Inverse Dynamics: */
603 void
604 StructDispNode::Update(const VectorHandler& X, InverseDynamics::Order iOrder)
605 {
606  integer iFirstIndex = iGetFirstIndex();
607  switch (iOrder) {
608  case InverseDynamics::POSITION:
609  XCurr = Vec3(X, iFirstIndex + 1);
610  break;
611 
612  case InverseDynamics::VELOCITY:
613  VCurr = Vec3(X, iFirstIndex + 1);
614  break;
615 
616  case InverseDynamics::ACCELERATION:
617  XPPCurr = Vec3(X, iFirstIndex + 1);
618  break;
619 
620  default:
621  ASSERT(0);
622  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
623  }
624 }
625 
626 /* Funzioni di inizializzazione, ereditate da DofOwnerOwner */
627 void
628 StructDispNode::SetInitialValue(VectorHandler& X)
629 {
630  /* FIXME: why is this called? */
631  integer iIndex = iGetFirstIndex();
632 
633  X.Put(iIndex + 1, XCurr);
634  X.Put(iIndex + 4, VCurr);
635 }
636 
637 
638 void
639 StructDispNode::SetValue(DataManager *pDM,
640  VectorHandler& X, VectorHandler& XP,
641  SimulationEntity::Hints *ph)
642 {
643 #ifdef MBDYN_X_RELATIVE_PREDICTION
644  // FIXME
645  if (pRefNode) {
646  Vec3 Xtmp = XPrev - pRefNode->GetXCurr();
647  const Mat3x3& R0 = pRefNode->GetRCurr();
648  const Vec3& V0 = pRefNode->GetVCurr();
649  const Vec3& W0 = pRefNode->GetWCurr();
650 
651  XPrev = R0.MulTV(Xtmp);
652  RPrev = R0.MulTM(RCurr);
653  VPrev = R0.MulTV(VCurr - V0 - W0.Cross(Xtmp));
654  WPrev = R0.MulTV(WCurr - W0);
655 
656 #if 0
657  std::cout << "StructNode(" << GetLabel() << "): "
658  "SetValue: X=" << XPrev
659  << ", R=" << RPrev
660  << ", V=" << VPrev
661  << ", W=" << WPrev
662  << std::endl;
663 #endif
664 
665  } else
666 #endif /* MBDYN_X_RELATIVE_PREDICTION */
667  {
668  /* FIXME: in any case, we start with Crank-Nicolson ... */
669  XPrev = XCurr;
670  VPrev = VCurr;
671  }
672 
673  integer iFirstIndex = iGetFirstIndex();
674  X.Put(iFirstIndex + 1, XPrev);
675  XP.Put(iFirstIndex + 1, VPrev);
676 }
677 
678 
679 void
680 StructDispNode::BeforePredict(VectorHandler& X,
681  VectorHandler& XP,
682  VectorHandler& XPr,
683  VectorHandler& XPPr) const
684 {
685 #ifdef MBDYN_X_RELATIVE_PREDICTION
686  // FIXME
687  integer iFirstPos = iGetFirstIndex();
688 
689  /* If pRefNode is defined, the prediction is made
690  * on the data in the reference frame it provides */
691  if (pRefNode) {
692 
693  /*
694  x_r = R_0^T * ( x - x_0 )
695  R_r = R_0^T * R
696  v_r = R_0^T * ( v - v_0 - omega_0 \times ( x - x_0 ) )
697  omega_r = R_0^T * ( omega - omega_0 )
698  */
699  Vec3 Xtmp = XCurr - pRefNode->GetXCurr();
700  const Mat3x3& R0 = pRefNode->GetRCurr();
701  const Vec3& V0 = pRefNode->GetVCurr();
702  const Vec3& W0 = pRefNode->GetWCurr();
703 
704  XCurr = R0.MulTV(Xtmp);
705  RCurr = R0.MulTM(RCurr);
706  VCurr = R0.MulTV(VCurr - V0 - W0.Cross(Xtmp));
707  WCurr = R0.MulTV(WCurr - W0);
708 
709  /* update state vectors with relative position and velocity */
710  X.Put(iFirstPos + 1, XCurr);
711  XP.Put(iFirstPos + 1, VCurr);
712  XPr.Put(iFirstPos + 1, XPrev);
713  XPPr.Put(iFirstPos + 1, VPrev);
714 
715 #if 0
716  std::cout << "StructNode(" << GetLabel() << "): "
717  "BeforePredict: X=" << XCurr
718  << ", R=" << RCurr
719  << ", V=" << VCurr
720  << ", W=" << WCurr
721  << std::endl;
722 #endif
723  }
724 #endif /* MBDYN_X_RELATIVE_PREDICTION */
725 
726  XPrev = XCurr;
727  VPrev = VCurr;
728 }
729 
730 void
731 StructDispNode::AfterPredict(VectorHandler& X, VectorHandler& XP)
732 {
733  integer iFirstIndex = iGetFirstIndex();
734 
735  /* Spostamento e velocita' aggiornati */
736  XCurr = Vec3(X, iFirstIndex + 1);
737  VCurr = Vec3(XP, iFirstIndex + 1);
738 
739 #ifdef MBDYN_X_RELATIVE_PREDICTION
740  if (pRefNode) {
741  // FIXME
742 
743  /*
744  x = x_0 + R_0 * x_r
745  R = R_0 * R_r
746  v = v_0 + omega_0 \times ( R_0 * x_r ) + R_0 * v_r
747  omega = omega_0 + R_0 * omega_r
748  */
749  Vec3 X0 = pRefNode->GetXCurr();
750  Mat3x3 R0 = pRefNode->GetRCurr();
751  Vec3 V0 = pRefNode->GetVCurr();
752  Vec3 W0 = pRefNode->GetWCurr();
753 
754  XCurr = R0*XCurr; /* temporary */
755  RCurr = R0*RCurr;
756  VCurr = V0 + W0.Cross(XCurr) + R0*VCurr;
757  WCurr = W0 + R0*WCurr;
758  XCurr += X0; /* plus reference */
759 
760  /* alcuni usano anche le predizioni dei parametri
761  * di rotazione e delle loro derivate come riferimento
762  * (approccio updated-updated); quindi calcolo
763  * i parametri di riferimento come i parametri
764  * che danno una predizione pari alla variazione
765  * di R0 piu' l'incremento relativo, e le derivate
766  * dei parametri corrispondenti */
767  gRef = Vec3(CGR_Rot::Param, R0*RDelta.MulMT(pRefNode->GetRPrev()));
768  gPRef = Mat3x3(CGR_Rot::MatGm1, gRef)*WCurr;
769 
770  /* to be safe, the correct values are put back
771  * in the state vectors */
772  X.Put(iFirstIndex + 1, XCurr);
773  XP.Put(iFirstIndex + 1, VCurr);
774 
775 #if 0
776  std::cout << "StructNode(" << GetLabel() << "): "
777  "AfterPredict: X=" << XCurr
778  << ", R=" << RCurr
779  << ", V=" << VCurr
780  << ", W=" << WCurr
781  << std::endl;
782 #endif
783  }
784 #endif /* MBDYN_X_RELATIVE_PREDICTION */
785 }
786 
787 /* Inverse Dynamics: */
788 void
789 StructDispNode::AfterConvergence(const VectorHandler& X,
790  const VectorHandler& XP,
791  const VectorHandler& XPP)
792 {
793  XPrev = XCurr;
794  VPrev = VCurr;
795  XPPPrev = XPPCurr;
796 }
797 
798 /*
799  * Metodi per l'estrazione di dati "privati".
800  * Si suppone che l'estrattore li sappia interpretare.
801  * Come default non ci sono dati privati estraibili
802  */
803 unsigned int
804 StructDispNode::iGetNumPrivData(void) const
805 {
806  unsigned i =
807  3 // X
808  + 3 // x (R^T * X)
809  + 3 // Phi
810  + 3 // XP
811  + 3 // xP (R^T * XP)
812  + 3 // Omega
813  + 3 // omega (R^T * Omega)
814  + 3 // Euler angles (123)
815  + 3 // Euler angles (313)
816  + 3 // Euler angles (321)
817  + 4; // Euler parameters
818 
819  if (bComputeAccelerations()) {
820  i +=
821  3 // XPP
822  + 3 // xPP (R^T * XPP)
823  + 3 // OmegaP
824  + 3; // omegaP (R^T * OmegaP)
825  }
826 
827  return i;
828 }
829 
830 /*
831  * Maps a string (possibly with substrings) to a private data;
832  * returns a valid index ( > 0 && <= iGetNumPrivData()) or 0
833  * in case of unrecognized data; error must be handled by caller
834  */
835 unsigned int
836 StructDispNode::iGetPrivDataIdx(const char *s) const
837 {
838  long idx;
839  char *next;
840  std::string sDataName(s);
841 
842  const char *brk = std::strchr(s, '[' /*]*/ );
843  if (brk == 0) {
844  return 0;
845  }
846 
847  size_t len = brk - s;;
848  brk++;
849 
850  errno = 0;
851  idx = strtol(brk, &next, 10);
852  int save_errno = errno;
853  if (next == brk || strcmp(next, /*[*/ "]") != 0) {
854  return 0;
855  }
856 
857  if (save_errno == ERANGE) {
858  silent_cerr("StructNode(" << GetLabel() << "): "
859  "warning, private data index "
860  << std::string(brk, next - brk)
861  << " overflows" << std::endl);
862  return 0;
863  }
864 
865  /*
866  X 0 + idx idx = {1,3}
867  x 3 + idx idx = {1,3}
868  Phi 6 + idx idx = {1,3}
869  XP 9 + idx idx = {1,3}
870  x 12 + idx idx = {1,3}
871  Omega 15 + idx idx = {1,3}
872  omega 18 + idx idx = {1,3}
873  E | E123 21 + idx idx = {1,3}
874  E313 24 + idx idx = {1,3}
875  E321 27 + idx idx = {1,3}
876  PE 31 + idx idx = {0,3}
877  -------------------------------------------
878  XPP 34 + idx idx = {1,3}
879  xPP 37 + idx idx = {1,3}
880  OmegaP 40 + idx idx = {1,3}
881  omegaP 43 + idx idx = {1,3}
882  */
883 
884  if (strncmp(s, "PE", len) == 0) {
885  if (idx < 0 || idx > 3) {
886  return 0;
887  }
888 
889  return 31 + idx;
890  }
891 
892  if (idx < 1 || idx > 3) {
893  return 0;
894  }
895 
896  if (strncmp(s, "X", len) == 0) {
897  return 0 + idx;
898  }
899 
900  if (strncmp(s, "x", len) == 0) {
901  return 3 + idx;
902  }
903 
904  if (strncmp(s, "Phi", len) == 0) {
905  return 6 + idx;
906  }
907 
908  if (strncmp(s, "XP", len) == 0) {
909  return 9 + idx;
910  }
911 
912  if (strncmp(s, "xP", len) == 0) {
913  return 12 + idx;
914  }
915 
916  if (strncmp(s, "Omega", len) == 0) {
917  return 15 + idx;
918  }
919 
920  if (strncmp(s, "omega", len) == 0) {
921  return 18 + idx;
922  }
923 
924  if (strncmp(s, "E", len) == 0
925  || strncmp(s, "E123", len) == 0)
926  {
927  return 21 + idx;
928  }
929 
930  if (strncmp(s, "E313", len) == 0) {
931  return 24 + idx;
932  }
933 
934  if (strncmp(s, "E321", len) == 0) {
935  return 27 + idx;
936  }
937 
938  bool bca = false;
939  unsigned i;
940  if (strncmp(s, "XPP", len) == 0) {
941  bca = true;
942  i = 34 + idx;
943 
944  } else if (strncmp(s, "xPP", len) == 0) {
945  bca = true;
946  i = 37 + idx;
947 
948  } else if (strncmp(s, "OmegaP", len) == 0) {
949  bca = true;
950  i = 40 + idx;
951 
952  } else if (strncmp(s, "omegaP", len) == 0) {
953  bca = true;
954  i = 43 + idx;
955 
956  } else {
957  // error
958  return 0;
959  }
960 
961  // NOTE: bComputeAccels is set only if iGetPrivDataIdx() is called
962  // first; it is not when the (deprecated) idx is directly used.
963  if (bca) {
964  if (!const_cast<StructDispNode *>(this)->ComputeAccelerations(true)) {
965  silent_cerr("StructNode(" << GetLabel() << "): "
966  "request to compute accelerations failed, requested by private data \"" << sDataName << "\""
967  << std::endl);
968  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
969  }
970  }
971 
972  return i;
973 }
974 
975 /*
976  * Returns the current value of a private data
977  * with 0 < i <= iGetNumPrivData()
978  */
979 doublereal
980 StructDispNode::dGetPrivData(unsigned int i) const
981 {
982  switch (i) {
983  case 1:
984  case 2:
985  case 3:
986  return XCurr(i);
987 
988  case 4:
989  case 5:
990  case 6:
991  return XCurr(i - 3);
992 
993  case 7:
994  case 8:
995  case 9: {
996  return 0.;
997  }
998 
999  case 10:
1000  case 11:
1001  case 12:
1002  return VCurr(i - 9);
1003 
1004  case 13:
1005  case 14:
1006  case 15:
1007  return VCurr(i - 12);
1008 
1009  case 16:
1010  case 17:
1011  case 18:
1012  return 0.;
1013 
1014  case 19:
1015  case 20:
1016  case 21:
1017  return 0.;
1018 
1019  case 22:
1020  case 23:
1021  case 24:
1022  return 0.;
1023 
1024  case 25:
1025  case 26:
1026  case 27:
1027  return 0.;
1028 
1029  case 28:
1030  case 29:
1031  case 30:
1032  return 0.;
1033 
1034  case 31:
1035  case 32:
1036  case 33:
1037  case 34:
1038  return 0.;
1039 
1040  case 35:
1041  case 36:
1042  case 37:
1043  ASSERT(bComputeAccelerations() == true);
1044  return XPPCurr(i - 34);
1045 
1046  case 38:
1047  case 39:
1048  case 40:
1049  ASSERT(bComputeAccelerations() == true);
1050  return XPPCurr(i - 37);
1051 
1052  case 41:
1053  case 42:
1054  case 43:
1055  ASSERT(bComputeAccelerations() == true);
1056  return 0.;
1057 
1058  case 44:
1059  case 45:
1060  case 46:
1061  ASSERT(bComputeAccelerations() == true);
1062  return 0.;
1063  }
1064 
1065  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1066 }
1067 
1068 /* StructNode - end */
1069 
1070 
1071 
1072 /* StructDispNode - end */
1073 
1074 /* DynamicStructDispNode - begin */
1075 
1076 DynamicStructDispNode::DynamicStructDispNode(unsigned int uL,
1077  const DofOwner* pDO,
1078  const Vec3& X0,
1079  const Vec3& V0,
1080  const StructNode *pRN,
1081  const RigidBodyKinematics *pRBK,
1082  doublereal dPosStiff,
1083  doublereal dVelStiff,
1084  OrientationDescription ood,
1085  flag fOut)
1086 :
1087 StructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut),
1088 bComputeAccels((fOut & OUTPUT_ACCELERATIONS) == OUTPUT_ACCELERATIONS),
1089 pAutoStr(0)
1090 {
1091  bOutputAccels = bComputeAccels;
1092 }
1093 
1094 /* Distruttore (per ora e' banale) */
1095 DynamicStructDispNode::~DynamicStructDispNode(void)
1096 {
1097  NO_OP;
1098 }
1099 
1100 StructDispNode::Type
1101 DynamicStructDispNode::GetStructDispNodeType(void) const
1102 {
1103  return StructDispNode::DYNAMIC;
1104 }
1105 
1106 /* rigid-body kinematics */
1107 const Vec3&
1108 DynamicStructDispNode::GetXPP(void) const
1109 {
1110  return GetXPPCurr();
1111 }
1112 
1113 std::ostream&
1114 DynamicStructDispNode::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
1115 {
1116  integer iIndex = iGetFirstIndex();
1117 
1118  StructDispNode::DescribeDof(out, prefix, bInitial);
1119 
1120  if (bInitial == false) {
1121  out
1122  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
1123  "momentum [Bx,By,Bz]" << std::endl;
1124  }
1125 
1126  return out;
1127 }
1128 
1129 void
1130 DynamicStructDispNode::DescribeDof(std::vector<std::string>& desc, bool bInitial, int i) const
1131 {
1132  if (bInitial || i == -1 || (i >= 0 && i < 3)) {
1133  StructDispNode::DescribeDof(desc, bInitial, i);
1134 
1135  if (bInitial || (i >= 0 && i < 3)) {
1136  return;
1137  }
1138  }
1139 
1140  if (i == -1) {
1141  desc.resize(6);
1142 
1143  } else {
1144  desc.resize(1);
1145  }
1146 
1147  std::ostringstream os;
1148  os << "StructDispNode(" << GetLabel() << ")";
1149 
1150  if (i == -1) {
1151  std::string name = os.str();
1152 
1153  for (i = 3; i < 6; i++) {
1154  os.str(name);
1155  os.seekp(0, std::ios_base::end);
1156  os << ": " << sdn_dof[i/3] << xyz[i%3];
1157  desc[i] = os.str();
1158  }
1159 
1160  } else {
1161  if (i < 3 || i >= 6) {
1162  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1163  }
1164 
1165  os << ": " << sdn_dof[i/3] << xyz[i%3];
1166  desc[0] = os.str();
1167  }
1168 }
1169 
1170 std::ostream&
1171 DynamicStructDispNode::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
1172 {
1173  if (bInitial == false) {
1174  integer iIndex = iGetFirstIndex();
1175 
1176  out
1177  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
1178  "momentum definition [Bx,By,Bz]" << std::endl;
1179  }
1180 
1181  StructDispNode::DescribeEq(out, prefix, bInitial);
1182 
1183  return out;
1184 }
1185 
1186 void
1187 DynamicStructDispNode::DescribeEq(std::vector<std::string>& desc, bool bInitial, int i) const
1188 {
1189  if (bInitial || i == -1 || (i >= 3 && i < 6)) {
1190  int new_i = i;
1191  if (!bInitial && i != -1) {
1192  new_i = i - 3;
1193  }
1194  StructDispNode::DescribeEq(desc, bInitial, new_i);
1195 
1196  if (bInitial || (i >= 3 && i < 6)) {
1197  return;
1198  }
1199  }
1200 
1201  if (i == -1) {
1202  desc.resize(6);
1203  for (int j = 0; j < 3; j++) {
1204  desc[3 + j] = desc[j];
1205  }
1206 
1207  } else {
1208  desc.resize(1);
1209  }
1210 
1211  std::ostringstream os;
1212  os << "StructDispNode(" << GetLabel() << ")";
1213 
1214  if (i == -1) {
1215  std::string name(os.str());
1216 
1217  for (i = 0; i < 3; i++) {
1218  os.str(name);
1219  os.seekp(0, std::ios_base::end);
1220  os << ": " << sdn_eq[i/3] << xyz[i%3];
1221  desc[i] = os.str();
1222  }
1223 
1224  } else {
1225  if (i < 0 || i >= 3) {
1226  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1227  }
1228 
1229  os << ": " << sdn_eq[i/3] << xyz[i%3];
1230  desc[0] = os.str();
1231  }
1232 }
1233 
1234 /* Usato dalle forze astratte, dai bulk ecc., per assemblare le forze
1235  * al posto giusto */
1236 integer
1237 DynamicStructDispNode::iGetFirstRowIndex(void) const
1238 {
1239  return iGetFirstMomentumIndex();
1240 }
1241 
1242 /* delegate to autostr node */
1243 void
1244 DynamicStructDispNode::AddInertia(const doublereal& dm) const
1245 {
1246  /* FIXME: do it only if to be output... */
1247  if (bComputeAccelerations()) {
1248  pAutoStr->AddInertia(dm);
1249  }
1250 }
1251 
1252 const Vec3&
1253 DynamicStructDispNode::GetBCurr(void) const
1254 {
1255  return pAutoStr->GetBCurr();
1256 }
1257 
1258 const Vec3&
1259 DynamicStructDispNode::GetBPCurr(void) const
1260 {
1261  return pAutoStr->GetBPCurr();
1262 }
1263 
1264 void
1265 DynamicStructDispNode::Update(const VectorHandler& X, const VectorHandler& XP)
1266 {
1267  StructDispNode::Update(X, XP);
1268  if (bComputeAccelerations()) {
1269  ASSERT(pAutoStr != 0);
1270 
1271  // FIXME: based on values set during previous
1272  // of AutomaticStructural::AssRes()
1273  pAutoStr->ComputeAccelerations(XPPCurr);
1274  }
1275 }
1276 
1277 void
1278 DynamicStructDispNode::AfterConvergence(const VectorHandler& X,
1279  const VectorHandler& XP)
1280 {
1281  if (bComputeAccelerations()) {
1282  ASSERT(pAutoStr != 0);
1283  pAutoStr->ComputeAccelerations(XPPCurr);
1284  }
1285 }
1286 
1287 void
1288 DynamicStructDispNode::BeforePredict(VectorHandler& X,
1289  VectorHandler& XP,
1290  VectorHandler& XPr,
1291  VectorHandler& XPPr) const
1292 {
1293  if (bComputeAccelerations()) {
1294  XPPPrev = XPPCurr;
1295  }
1296 
1297  StructDispNode::BeforePredict(X, XP, XPr, XPPr);
1298 }
1299 
1300 /* Restituisce il valore del dof iDof;
1301  * se differenziale, iOrder puo' essere = 1 per la derivata */
1302 const doublereal&
1303 DynamicStructDispNode::dGetDofValue(int iDof, int iOrder) const
1304 {
1305  ASSERT(iDof >= 1 && iDof <= 3);
1306  ASSERT(iOrder >= 0 && iOrder <= 2);
1307 
1308  if (iOrder == 2) {
1309  /* FIXME: should not happen */
1310  ASSERT(bComputeAccelerations());
1311  if (!bComputeAccelerations()) {
1312  silent_cerr("DynamicStructDispNode::dGetDofValue("
1313  << iDof << "," << iOrder << "): "
1314  "accelerations are not computed while they should"
1315  << std::endl);
1316  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1317  }
1318 
1319 #if 1
1320  /* FIXME: might need to compute them in order to be
1321  * as up to date as possible; however, elements that contribute
1322  * to inertia should assemble first...
1323  */
1324  pAutoStr->ComputeAccelerations(XPPCurr);
1325 #endif
1326 
1327  return XPPCurr(iDof);
1328  }
1329 
1330  return StructDispNode::dGetDofValue(iDof, iOrder);
1331 }
1332 
1333 /* Restituisce il valore del dof iDof al passo precedente;
1334  * se differenziale, iOrder puo' essere = 1 per la derivata */
1335 const doublereal&
1336 DynamicStructDispNode::dGetDofValuePrev(int iDof, int iOrder) const
1337 {
1338  ASSERT(iDof >= 1 && iDof <= 3);
1339  ASSERT(iOrder == 0 || iOrder == 1);
1340 
1341  if (iOrder == 2) {
1342  /* FIXME: should not happen */
1343  ASSERT(bComputeAccelerations());
1344  if (!bComputeAccelerations()) {
1345  silent_cerr("DynamicStructDispNode::dGetDofValuePrev("
1346  << iDof << "," << iOrder << "): "
1347  "accelerations are not computed while they should"
1348  << std::endl);
1349  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1350  }
1351 
1352  return XPPPrev(iDof);
1353  }
1354 
1355  return StructDispNode::dGetDofValuePrev(iDof, iOrder);
1356 }
1357 
1358 /* Setta il valore del dof iDof a dValue;
1359  * se differenziale, iOrder puo' essere = 1 per la derivata */
1360 void
1361 DynamicStructDispNode::SetDofValue(const doublereal& dValue,
1362  unsigned int iDof,
1363  unsigned int iOrder /* = 0 */ )
1364 {
1365  ASSERT(iDof >= 1 && iDof <= 3);
1366  ASSERT(iOrder == 0 || iOrder == 1);
1367 
1368  if (iOrder == 2) {
1369  /* FIXME: should not happen */
1370  ASSERT(bComputeAccelerations());
1371  if (!bComputeAccelerations()) {
1372  silent_cerr("DynamicStructNode::SetDofValue("
1373  << dValue << "," << iDof << "," << iOrder << "): "
1374  "accelerations are not computed while they should"
1375  << std::endl);
1376  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1377  }
1378 
1379  XPPCurr(iDof) = dValue;
1380 
1381  } else {
1382  StructDispNode::SetDofValue(iDof, iOrder);
1383  }
1384 }
1385 
1386 bool
1387 DynamicStructDispNode::ComputeAccelerations(bool b)
1388 {
1389  bComputeAccels = b;
1390  return true;
1391 }
1392 
1393 void
1394 DynamicStructDispNode::SetOutputFlag(flag f)
1395 {
1396  if (f & StructDispNode::OUTPUT_ACCELERATIONS) {
1397  // ignore result
1398  ComputeAccelerations(true);
1399  }
1400  ToBeOutput::SetOutputFlag(f);
1401 }
1402 
1403 /* DynamicStructDispNode - end */
1404 
1405 /* StaticStructDispNode - begin */
1406 
1407 StaticStructDispNode::StaticStructDispNode(unsigned int uL,
1408  const DofOwner* pDO,
1409  const Vec3& X0,
1410  const Vec3& V0,
1411  const StructNode *pRN,
1412  const RigidBodyKinematics *pRBK,
1413  doublereal dPosStiff,
1414  doublereal dVelStiff,
1415  OrientationDescription ood,
1416  flag fOut)
1417 :
1418 StructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut)
1419 {
1420  NO_OP;
1421 }
1422 
1423 StaticStructDispNode::~StaticStructDispNode(void)
1424 {
1425  NO_OP;
1426 }
1427 
1428 StructDispNode::Type
1429 StaticStructDispNode::GetStructDispNodeType(void) const
1430 {
1431  return StructDispNode::STATIC;
1432 }
1433 
1434 /* StaticStructDispNode - end */
1435 
1436 /* StructNode - begin */
1437 
1438 /* Costruttore definitivo */
1439 StructNode::StructNode(unsigned int uL,
1440  const DofOwner* pDO,
1441  const Vec3& X0,
1442  const Mat3x3& R0,
1443  const Vec3& V0,
1444  const Vec3& W0,
1445  const StructNode *pRN,
1446  const RigidBodyKinematics *pRBK,
1447  doublereal dPosStiff,
1448  doublereal dVelStiff,
1449  bool bOmRot,
1450  OrientationDescription ood,
1451  flag fOut)
1452 : StructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut),
1453 RPrev(R0),
1454 RRef(R0),
1455 RCurr(R0),
1456 gRef(Zero3),
1457 gCurr(Zero3),
1458 gPRef(Zero3),
1459 gPCurr(Zero3),
1460 WPrev(W0),
1461 WRef(W0),
1462 WCurr(W0),
1463 WPCurr(Zero3),
1464 WPPrev(Zero3),
1465 bOmegaRot(bOmRot)
1466 #ifdef USE_AUTODIFF
1467 ,bUpdateRotation(true)
1468 ,dCoefGrad(0.) // This should be safe because the time step should never be zero
1469 #endif
1470 {
1471  NO_OP;
1472 }
1473 
1474 /* Distruttore (per ora e' banale) */
1475 StructNode::~StructNode(void)
1476 {
1477  NO_OP;
1478 }
1479 
1480 const Mat3x3&
1481 StructNode::GetR(void) const
1482 {
1483  return GetRCurr();
1484 }
1485 
1486 const Vec3&
1487 StructNode::GetW(void) const
1488 {
1489  return GetWCurr();
1490 }
1491 
1492 const Vec3&
1493 StructNode::GetWP(void) const
1494 {
1495  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1496 }
1497 
1498 std::ostream&
1499 StructNode::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
1500 {
1501  integer iIndex = iGetFirstIndex();
1502 
1503  out
1504  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
1505  "position [px,py,pz]" << std::endl
1506  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
1507  "orientation parameters [gx,gy,gz]" << std::endl;
1508 
1509  if (bInitial) {
1510  iIndex += 6;
1511  out
1512  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
1513  "linear velocity [vx,vy,vz]" << std::endl
1514  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
1515  "angular velocity [wx,wy,wz]" << std::endl;
1516  }
1517 
1518  return out;
1519 }
1520 
1521 void
1522 StructNode::DescribeDof(std::vector<std::string>& desc, bool bInitial, int i) const
1523 {
1524  if (i == -1) {
1525  if (bInitial) {
1526  desc.resize(12);
1527 
1528  } else {
1529  desc.resize(6);
1530  }
1531 
1532  } else {
1533  desc.resize(1);
1534  }
1535 
1536  std::ostringstream os;
1537  os << "StructNode(" << GetLabel() << ")";
1538 
1539  // always uses initial_dof[] becuase dof[]
1540  // and initial_dof[] are the same up to 6
1541  int iend = bInitial ? 12 : 6;
1542  if (i == -1) {
1543  std::string name = os.str();
1544 
1545  for (i = 0; i < iend; i++) {
1546  os.str(name);
1547  os.seekp(0, std::ios_base::end);
1548  os << ": " << sn_initial_dof[i/3] << xyz[i%3];
1549  desc[i] = os.str();
1550  }
1551 
1552  } else {
1553  if (i < 0 || i >= iend) {
1554  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1555  }
1556  os << ": " << sn_initial_dof[i/3] << xyz[i%3];
1557  desc[0] = os.str();
1558  }
1559 }
1560 
1561 std::ostream&
1562 StructNode::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
1563 {
1564  integer iIndex = iGetFirstIndex();
1565 
1566  if (bInitial) {
1567  out
1568  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
1569  "position [Px,Py,Pz]" << std::endl
1570  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
1571  "orientation [gx,gy,gz]" << std::endl
1572  << prefix << iIndex + 7 << "->" << iIndex + 9 << ": "
1573  "linear velocity [vx,vy,vz]" << std::endl
1574  << prefix << iIndex + 10 << "->" << iIndex + 12 << ": "
1575  "angular velocity [wx,wy,wz]" << std::endl;
1576  } else {
1577  if (dynamic_cast<const DynamicStructNode*>(this) != 0
1578  || dynamic_cast<const ModalNode*>(this) != 0)
1579  {
1580  iIndex += 6;
1581  }
1582 
1583  out
1584  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
1585  "force equilibrium [Fx,Fy,Fz]" << std::endl
1586  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
1587  "moment equilibrium [Mx,My,Mz]" << std::endl;
1588  }
1589 
1590  return out;
1591 }
1592 
1593 void
1594 StructNode::DescribeEq(std::vector<std::string>& desc, bool bInitial, int i) const
1595 {
1596  if (i == -1) {
1597  if (bInitial) {
1598  desc.resize(12);
1599 
1600  } else {
1601  desc.resize(6);
1602  }
1603 
1604  } else {
1605  desc.resize(1);
1606  }
1607 
1608  std::ostringstream os;
1609  os << "StructNode(" << GetLabel() << ")";
1610 
1611  if (i == -1) {
1612  std::string name(os.str());
1613 
1614  if (bInitial) {
1615  for (i = 0; i < 12; i++) {
1616  os.str(name);
1617  os.seekp(0, std::ios_base::end);
1618  os << ": " << sn_initial_eq[i/3] << xyz[i%3];
1619  desc[i] = os.str();
1620  }
1621 
1622  } else {
1623  for (i = 0; i < 6; i++) {
1624  os.str(name);
1625  os.seekp(0, std::ios_base::end);
1626  os << ": " << sn_eq[2 + i/3] << xyz[i%3];
1627  desc[i] = os.str();
1628  }
1629  }
1630 
1631  } else {
1632  if (bInitial) {
1633  if (i < 0 || i >= 12) {
1634  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1635  }
1636 
1637  os << ": " << sn_initial_eq[i/3] << xyz[i%3];
1638 
1639  } else {
1640  if (i < 0 || i >= 6) {
1641  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1642  }
1643 
1644  os << ": " << sn_eq[2 + i/3] << xyz[i%3];
1645  }
1646  desc[0] = os.str();
1647  }
1648 }
1649 
1650 /* Contributo del nodo strutturale al file di restart */
1651 std::ostream&
1652 StructNode::Restart(std::ostream& out) const
1653 {
1654  out << " structural: " << GetLabel() << ", ";
1655  if (GetStructNodeType() == StructNode::DYNAMIC) {
1656  out << "dynamic";
1657  } else if (GetStructNodeType() == StructNode::STATIC) {
1658  out << "static";
1659  }
1660  out << ", reference, global, ";
1661  XCurr.Write(out, ", ")
1662  << ", reference, global, 1, ", (RCurr.GetVec(1)).Write(out, ", ")
1663  << ", 2, ", (RCurr.GetVec(2)).Write(out, ", ")
1664  << ", reference, global, ",
1665  VCurr.Write(out, ", ")
1666  << ", reference, global, ",
1667  WCurr.Write(out, ", ") << ", assembly, "
1668  << dPositionStiffness << ", "
1669  << dVelocityStiffness << ", "
1670  << bOmegaRot
1671  << ", scale, " << pGetDofOwner()->dGetScale() << ';' << std::endl;
1672 
1673  return out;
1674 }
1675 
1676 
1677 /* Restituisce il valore del dof iDof;
1678  * se differenziale, iOrder puo' essere = 1 per la derivata */
1679 const doublereal&
1680 StructNode::dGetDofValue(int iDof, int iOrder) const
1681 {
1682  ASSERT(iDof >= 1 && iDof <= 6);
1683  ASSERT(iOrder == 0 || iOrder == 1);
1684  if (iDof >= 1 && iDof <= 3) {
1685  if (iOrder == 0) {
1686  return XCurr(iDof);
1687  } else if (iOrder == 1) {
1688  return VCurr(iDof);
1689  }
1690  } else if (iDof >= 4 && iDof <= 6) {
1691  if (iOrder == 1) {
1692  return WCurr(iDof - 3);
1693  } else if (iOrder == 0) {
1694  silent_cerr("StructNode(" << GetLabel() << "): "
1695  "unable to return angles" << std::endl);
1696  throw StructNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
1697  }
1698  } else {
1699  silent_cerr("StructNode(" << GetLabel() << "): "
1700  "required dof " << iDof << " (order " << iOrder << ") "
1701  "is not available." << std::endl);
1702  throw StructNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
1703  }
1704 
1705  /* dummy return value to workaround compiler complains */
1706  static doublereal dmy = 0.;
1707  return dmy;
1708 }
1709 
1710 /* Restituisce il valore del dof iDof al passo precedente;
1711  * se differenziale, iOrder puo' essere = 1 per la derivata */
1712 const doublereal&
1713 StructNode::dGetDofValuePrev(int iDof, int iOrder) const
1714 {
1715  ASSERT(iDof >= 1 && iDof <= 6);
1716  ASSERT(iOrder == 0 || iOrder == 1);
1717  if (iDof >= 1 && iDof <= 3) {
1718  if (iOrder == 0) {
1719  return XPrev(iDof);
1720  } else if (iOrder == 1) {
1721  return VPrev(iDof);
1722  }
1723  } else if (iDof >= 4 && iDof <= 6) {
1724  if (iOrder == 1) {
1725  return WPrev(iDof - 3);
1726  } else if (iOrder == 0) {
1727  silent_cerr("StructNode(" << GetLabel() << "): "
1728  "unable to return angles" << std::endl);
1729  throw StructNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
1730  }
1731  } else {
1732  silent_cerr("StructNode(" << GetLabel() << "): "
1733  "required dof " << iDof << " (order " << iOrder << ") "
1734  "is not available." << std::endl);
1735  throw StructNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
1736  }
1737 
1738  /* dummy return value to workaround compiler complains */
1739  static doublereal dmy = 0.;
1740  return dmy;
1741 }
1742 
1743 /* Setta il valore del dof iDof a dValue;
1744  * se differenziale, iOrder puo' essere = 1 per la derivata */
1745 void
1746 StructNode::SetDofValue(const doublereal& dValue,
1747  unsigned int iDof,
1748  unsigned int iOrder /* = 0 */ )
1749 {
1750  ASSERT(iDof >= 1 && iDof <= 6);
1751  ASSERT(iOrder == 0 || iOrder == 1);
1752  if (iDof >= 1 && iDof <= 3) {
1753  if (iOrder == 0) {
1754  XCurr(iDof) = dValue;
1755 
1756  } else if (iOrder == 1) {
1757  VCurr(iDof) = dValue;
1758  }
1759 
1760  } else if (iDof >= 4 && iDof <= 6) {
1761  if (iOrder == 1) {
1762  WCurr(iDof - 3) = dValue;
1763 
1764  } else if (iOrder == 0) {
1765  silent_cerr("StructNode(" << GetLabel() << "): "
1766  "unable to set angles" << std::endl);
1767  throw StructNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
1768  }
1769 
1770  } else {
1771  silent_cerr("StructNode(" << GetLabel() << "): "
1772  "required dof " << iDof << " (order " << iOrder << ") "
1773  "is not available." << std::endl);
1774  throw StructNode::ErrGeneric(MBDYN_EXCEPT_ARGS);
1775  }
1776 }
1777 
1778 void
1779 StructNode::OutputPrepare(OutputHandler &OH)
1780 {
1781  if (bToBeOutput()) {
1782 #ifdef USE_NETCDF
1783  if (OH.UseNetCDF(OutputHandler::STRNODES)) {
1784  ASSERT(OH.IsOpen(OutputHandler::NETCDF));
1785 
1786  // node
1787  const char *type;
1788  switch (GetStructNodeType()) {
1789  case STATIC:
1790  type = "static";
1791  break;
1792 
1793  case DYNAMIC:
1794  type = "dynamic";
1795  break;
1796 
1797  case MODAL:
1798  type = "modal";
1799  break;
1800 
1801  case DUMMY:
1802  type = "dummy";
1803  break;
1804 
1805  default:
1806  pedantic_cerr("StructNode::OutputPrepare(" << GetLabel() << "): "
1807  "warning, unknown node type?" << std::endl);
1808  type = "unknown";
1809  break;
1810  }
1811 
1812  std::ostringstream os;
1813  os << "node.struct." << GetLabel();
1814  (void)OH.CreateVar(os.str(), type);
1815 
1816  // node sub-data
1817  os << '.';
1818  std::string name = os.str();
1819 
1820  Var_X = OH.CreateVar<Vec3>(name + "X", "m",
1821  "global position vector (X, Y, Z)");
1822 
1823  Var_Phi = OH.CreateRotationVar(name, "", od, "global");
1824 
1825  Var_XP = OH.CreateVar<Vec3>(name + "XP", "m/s",
1826  "global velocity vector (v_X, v_Y, v_Z)");
1827 
1828  Var_Omega = OH.CreateVar<Vec3>(name + "Omega", "radian/s",
1829  "global angular velocity vector (omega_X, omega_Y, omega_Z)");
1830 
1831  // accelerations
1832  if (bOutputAccels) {
1833  Var_XPP = OH.CreateVar<Vec3>(name + "XPP", "m/s^2",
1834  "global acceleration vector (a_X, a_Y, a_Z)");
1835 
1836  Var_OmegaP = OH.CreateVar<Vec3>(name + "OmegaP", "radian/s^2",
1837  "global angular acceleration vector (omegaP_X, omegaP_Y, omegaP_Z)");
1838  }
1839  }
1840 #endif // USE_NETCDF
1841  }
1842 }
1843 
1844 /* Output del nodo strutturale (da mettere a punto) */
1845 void
1846 StructNode::Output(OutputHandler& OH) const
1847 {
1848  if (bToBeOutput()) {
1849  Vec3 E;
1850  switch (od) {
1851  case EULER_123:
1852  E = MatR2EulerAngles123(RCurr)*dRaDegr;
1853  break;
1854 
1855  case EULER_313:
1856  E = MatR2EulerAngles313(RCurr)*dRaDegr;
1857  break;
1858 
1859  case EULER_321:
1860  E = MatR2EulerAngles321(RCurr)*dRaDegr;
1861  break;
1862 
1863  case ORIENTATION_VECTOR:
1864  E = RotManip::VecRot(RCurr);
1865  break;
1866 
1867  case ORIENTATION_MATRIX:
1868  break;
1869 
1870  default:
1871  /* impossible */
1872  break;
1873  }
1874 
1875 #ifdef USE_NETCDF
1876  if (OH.UseNetCDF(OutputHandler::STRNODES)) {
1877  Var_X->put_rec(XCurr.pGetVec(), OH.GetCurrentStep());
1878  switch (od) {
1879  case EULER_123:
1880  case EULER_313:
1881  case EULER_321:
1882  case ORIENTATION_VECTOR:
1883  Var_Phi->put_rec(E.pGetVec(), OH.GetCurrentStep());
1884  break;
1885 
1886  case ORIENTATION_MATRIX:
1887  Var_Phi->put_rec(RCurr.pGetMat(), OH.GetCurrentStep());
1888  break;
1889 
1890  default:
1891  /* impossible */
1892  break;
1893  }
1894  Var_XP->put_rec(VCurr.pGetVec(), OH.GetCurrentStep());
1895  Var_Omega->put_rec(WCurr.pGetVec(), OH.GetCurrentStep());
1896 
1897  if (bOutputAccels) {
1898  Var_XPP->put_rec(XPPCurr.pGetVec(), OH.GetCurrentStep());
1899  Var_OmegaP->put_rec(WPCurr.pGetVec(), OH.GetCurrentStep());
1900  }
1901  }
1902 #endif /* USE_NETCDF */
1903 
1904  if (OH.UseText(OutputHandler::STRNODES)) {
1905  std::ostream& out = OH.StrNodes();
1906  out
1907  << std::setw(8) << GetLabel()
1908  << " " << XCurr << " ";
1909  switch (od) {
1910  case EULER_123:
1911  case EULER_313:
1912  case EULER_321:
1913  case ORIENTATION_VECTOR:
1914  OH.StrNodes() << E;
1915  break;
1916 
1917  case ORIENTATION_MATRIX:
1918  OH.StrNodes() << RCurr;
1919  break;
1920 
1921  default:
1922  /* impossible */
1923  break;
1924  }
1925  OH.StrNodes() << " " << VCurr << " " << WCurr;
1926 
1927  if (bOutputAccels) {
1928  out
1929  << " " << XPPCurr
1930  << " " << WPCurr;
1931  }
1932  out << std::endl;
1933  }
1934  }
1935 }
1936 
1937 /* Aggiorna dati in base alla soluzione */
1938 void
1939 StructNode::Update(const VectorHandler& X, const VectorHandler& XP)
1940 {
1941 #ifdef USE_AUTODIFF
1942  bUpdateRotation = true;
1943 #endif
1944 
1945  integer iFirstIndex = iGetFirstIndex();
1946 
1947  XCurr = Vec3(X, iFirstIndex + 1);
1948  VCurr = Vec3(XP, iFirstIndex + 1);
1949 
1950  /* Nota: i g, gP non vengono incrementati */
1951  gCurr = Vec3(X, iFirstIndex + 4);
1952  gPCurr = Vec3(XP, iFirstIndex + 4);
1953 
1954 #if 0
1955  // test amplitude of orientation increment
1956  if (gCurr.Norm() > 1.) {
1957  silent_cerr("StructNode(" << GetLabel() << "): "
1958  "incremental rotation too large, YMMV" << std::endl);
1959  }
1960 #endif
1961 
1962  /* Matrice RDelta, incremento di rotazione da predetto a corrente;
1963  * Questo e' piu' efficiente */
1964  Mat3x3 RDelta(CGR_Rot::MatR, gCurr);
1965 
1966 #if 0
1967  /* Questo e' meno efficiente anche se sembra piu' elegante.
1968  * Il problema e' che per scrivere il manipolatore in forma
1969  * elegante bisogna aggiungere alla matrice le informazioni
1970  * di memorizzazione della funzione di manipolazione.
1971  * Oppure occorre un operatore ternario */
1972  RDelta = MatR << gCurr;
1973 #endif
1974 
1975  /* La matrice di rotazione corrente e' data dalla matrice predetta
1976  * (costante) moltiplicata per l'incremento totale occorso;
1977  * la velocita' angolare e' data dalla parte incrementale totale
1978  * piu' il contributo della velocita' di riferimento (costante) */
1979  RCurr = RDelta*RRef;
1980  WCurr = Mat3x3(CGR_Rot::MatG, gCurr)*gPCurr + RDelta*WRef;
1981 
1982 #if 0
1983  /* Nuovo manipolatore (forse e' meno efficiente) */
1984  WCurr = (CGR_Rot::MatG << gCurr)*gPCurr+RDelta*WRef;
1985 #endif
1986 }
1987 
1988 #ifdef USE_AUTODIFF
1989 inline void StructNode::GetgCurr(grad::Vector<grad::Gradient<iNumADVars>, 3>& g, doublereal dCoef, enum grad::FunctionCall func) const {
1990  using namespace grad;
1991 
1992  integer iFirstIndexLocal;
1993 
1994  switch (func) {
1995  case INITIAL_ASS_JAC:
1996  GRADIENT_ASSERT(dCoef == 1.);
1997  case INITIAL_DER_JAC:
1998  case REGULAR_JAC:
1999  iFirstIndexLocal = -1; // Always start from zero in order to save memory
2000  break;
2001 
2002  default:
2003  GRADIENT_ASSERT(false);
2004  }
2005 
2006  for (index_type i = 1; i <= 3; ++i) {
2007  Gradient<iNumADVars>& g_i = g(i);
2008  g_i.SetValuePreserve(gCurr(i));
2009  g_i.DerivativeResizeReset(0,
2010  iFirstIndexLocal + i,
2011  MapVectorBase::LOCAL,
2012  -dCoef);
2013  }
2014 }
2015 
2016 inline void StructNode::GetgPCurr(grad::Vector<grad::Gradient<iNumADVars>, 3>& gP, doublereal dCoef, enum grad::FunctionCall func) const {
2017  using namespace grad;
2018 
2019  integer iFirstIndexLocal;
2020 
2021  switch (func) {
2022  case INITIAL_ASS_JAC:
2023  GRADIENT_ASSERT(dCoef == 1.);
2024  case INITIAL_DER_JAC:
2025  case REGULAR_JAC:
2026  iFirstIndexLocal = -1; // Always start from zero in order to save memory
2027  break;
2028 
2029  default:
2030  GRADIENT_ASSERT(false);
2031  }
2032 
2033  // gPCurr is unused during initial assembly
2034  const Vec3& gPCurr = (func == INITIAL_ASS_JAC) ? WCurr : this->gPCurr;
2035 
2036  for (index_type i = 1; i <= 3; ++i) {
2037  Gradient<iNumADVars>& gP_i = gP(i);
2038  gP_i.SetValuePreserve(gPCurr(i));
2039  gP_i.DerivativeResizeReset(0,
2040  iFirstIndexLocal + i,
2041  MapVectorBase::LOCAL,
2042  -1.);
2043  }
2044 }
2045 
2046 template <typename T>
2047 void StructNode::UpdateRotation(const Mat3x3& RRef, const Vec3& WRef, const grad::Vector<T, 3>& g, const grad::Vector<T, 3>& gP, grad::Matrix<T, 3, 3>& R, grad::Vector<T, 3>& W, enum grad::FunctionCall func) const
2048 {
2049  using namespace grad;
2050 
2051  const T d = 4. / (4. + Dot(g, g));
2052 
2053  Matrix<T, 3, 3> RDelta;
2054 
2055  const T tmp1 = -g(3) * g(3);
2056  const T tmp2 = -g(2) * g(2);
2057  const T tmp3 = -g(1) * g(1);
2058  const T tmp4 = g(1) * g(2) * 0.5;
2059  const T tmp5 = g(2) * g(3) * 0.5;
2060  const T tmp6 = g(1) * g(3) * 0.5;
2061 
2062  RDelta(1,1) = (tmp1 + tmp2) * d * 0.5 + 1;
2063  RDelta(1,2) = (tmp4 - g(3)) * d;
2064  RDelta(1,3) = (tmp6 + g(2)) * d;
2065  RDelta(2,1) = (g(3) + tmp4) * d;
2066  RDelta(2,2) = (tmp1 + tmp3) * d * 0.5 + 1.;
2067  RDelta(2,3) = (tmp5 - g(1)) * d;
2068  RDelta(3,1) = (tmp6 - g(2)) * d;
2069  RDelta(3,2) = (tmp5 + g(1)) * d;
2070  RDelta(3,3) = (tmp2 + tmp3) * d * 0.5 + 1.;
2071 
2072  R = RDelta * RRef;
2073 
2074  switch (func) {
2075  case INITIAL_ASS_JAC:
2076  W = gP; // Note gP must be equal to WCurr during the initial assembly phase
2077  break;
2078 
2079  case INITIAL_DER_JAC:
2080  case REGULAR_JAC:
2081  {
2082  Matrix<T, 3, 3> G;
2083 
2084  const T tmp7 = 0.5 * g(1) * d;
2085  const T tmp8 = 0.5 * g(2) * d;
2086  const T tmp9 = 0.5 * g(3) * d;
2087 
2088  G(1,1) = d;
2089  G(1,2) = -tmp9;
2090  G(1,3) = tmp8;
2091  G(2,1) = tmp9;
2092  G(2,2) = d;
2093  G(2,3) = -tmp7;
2094  G(3,1) = -tmp8;
2095  G(3,2) = tmp7;
2096  G(3,3) = d;
2097 
2098  W = G * gP + RDelta * WRef; // Note that the first index of gP and g must be the same in order to work!
2099  }
2100  break;
2101 
2102  default:
2103  GRADIENT_ASSERT(false);
2104  }
2105 
2106 }
2107 
2108 void StructNode::UpdateRotation(doublereal dCoef, enum grad::FunctionCall func) const
2109 {
2110  using namespace grad;
2111 
2112  if (bUpdateRotation || dCoef != dCoefGrad) {
2113 #ifdef USE_MULTITHREAD
2114  if (!gradInUse.bIsInUse()) {
2115  // Another thread has locked this node
2116  // Wait until it is finished
2117 
2118  while (!gradInUse.bIsInUse()) {
2119  DEBUGCOUT("StructNode::UpdateRotation: StructNode(" << GetLabel() << ") is in use!\n");
2120  }
2121 
2122  if (!bUpdateRotation && dCoef == dCoefGrad) {
2123  gradInUse.ReSetInUse();
2124  return;
2125  } else {
2126  ASSERT(0);
2127  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2128  }
2129  }
2130 
2131  try {
2132 #else
2133  {
2134 #endif
2135  dCoefGrad = dCoef;
2136  Vector<Gradient<iNumADVars>, 3> gCurr_grad, gPCurr_grad;
2137 
2138  GetgCurr(gCurr_grad, dCoef, func);
2139  GetgPCurr(gPCurr_grad, dCoef, func);
2140 
2141  UpdateRotation(RRef, WRef, gCurr_grad, gPCurr_grad, RCurr_grad, WCurr_grad, func);
2142 
2143 #if GRADIENT_DEBUG > 0
2144  {
2145  const double dTol = 10 * std::numeric_limits<doublereal>::epsilon();
2146 
2147  Mat3x3 RDelta(CGR_Rot::MatR, gCurr);
2148 
2149  Mat3x3 RCurr_tmp = RDelta * RRef;
2150 
2151  bool bErr = false;
2152 
2153  for (index_type i = 1; i <= 3; ++i) {
2154  for (index_type j = 1; j <= 3; ++j) {
2155  if (std::abs(RCurr_grad(i, j).dGetValue() - RCurr_tmp(i, j)) > dTol) {
2156  bErr = true;
2157  }
2158  }
2159  }
2160 
2161  for (index_type i = 1; i <= 3; ++i) {
2162  for (index_type j = 1; j <= 3; ++j) {
2163  if (std::abs(RCurr_grad(i, j).dGetValue() - RCurr(i, j)) > dTol) {
2164  bErr = true;
2165  }
2166  }
2167 
2168  if (std::abs(WCurr_grad(i).dGetValue() - WCurr(i)) > dTol) {
2169  bErr = true;
2170  }
2171  }
2172 
2173  if (bErr) {
2174  std::cerr << "gCurr=" << gCurr << std::endl;
2175  std::cerr << "RCurr=" << std::endl;
2176  for (integer i = 1; i <= 3; ++i) {
2177  for (integer j = 1; j <= 3; ++j) {
2178  std::cerr << RCurr(i, j) << " ";
2179  }
2180  std::cerr << std::endl;
2181  }
2182 
2183  std::cerr << "RCurr_grad=" << std::endl;
2184 
2185  std::cerr << "RRef=" << std::endl;
2186  for (integer i = 1; i <= 3; ++i) {
2187  for (integer j = 1; j <= 3; ++j) {
2188  std::cerr << RRef(i, j) << " ";
2189  }
2190  std::cerr << std::endl;
2191  }
2192 
2193  std::cerr << "RPrev=" << std::endl;
2194  for (integer i = 1; i <= 3; ++i) {
2195  for (integer j = 1; j <= 3; ++j) {
2196  std::cerr << RPrev(i, j) << " ";
2197  }
2198  std::cerr << std::endl;
2199  }
2200 
2201  std::cerr << "RCurr_grad=" << std::endl;
2202 
2203  for (integer i = 1; i <= 3; ++i) {
2204  for (integer j = 1; j <= 3; ++j) {
2205  std::cerr << RCurr_grad(i, j).dGetValue() << " ";
2206  }
2207  std::cerr << std::endl;
2208  }
2209 
2210  std::cerr << "WCurr=" << WCurr << std::endl;
2211  std::cerr << "WCurr_grad=";
2212  for (integer i = 1; i <= 3; ++i) {
2213  std::cerr << WCurr_grad(i).dGetValue() << " ";
2214  }
2215  std::cerr << std::endl;
2216  GRADIENT_ASSERT(false);
2217  }
2218  }
2219 #endif
2220 
2221  bUpdateRotation = false;
2222 
2223 #if USE_MULTITHREAD
2224  } catch (...) {
2225  // Make sure that the spin lock will be reset in order to avoid infinite loops
2226  gradInUse.ReSetInUse();
2227  throw;
2228  }
2229 
2230  gradInUse.ReSetInUse();
2231 #else
2232  }
2233 #endif
2234  }
2235 }
2236 #endif
2237 
2238 /* Aggiorna dati in base alla soluzione */
2239 void
2240 StructNode::DerivativesUpdate(const VectorHandler& X, const VectorHandler& XP)
2241 {
2242 #ifdef USE_AUTODIFF
2243  bUpdateRotation = true;
2244 #endif
2245 
2246  integer iFirstIndex = iGetFirstIndex();
2247 
2248  /* Forza configurazione e velocita' al valore iniziale */
2249  const_cast<VectorHandler &>(X).Put(iFirstIndex + 1, XCurr);
2250  const_cast<VectorHandler &>(X).Put(iFirstIndex + 4, gCurr);
2251  const_cast<VectorHandler &>(XP).Put(iFirstIndex + 1, VCurr);
2252  const_cast<VectorHandler &>(XP).Put(iFirstIndex + 4, gPCurr);
2253 }
2254 
2255 
2256 /* Aggiorna dati in base alla soluzione durante l'assemblaggio iniziale */
2257 void
2258 StructNode::InitialUpdate(const VectorHandler& X)
2259 {
2260 #ifdef USE_AUTODIFF
2261  bUpdateRotation = true;
2262 #endif
2263 
2264  integer iFirstIndex = iGetFirstIndex();
2265 
2266  XCurr = Vec3(X, iFirstIndex + 1);
2267  VCurr = Vec3(X, iFirstIndex + 7);
2268 
2269  /* Nota: g viene incrementato */
2270  gCurr = Vec3(X, iFirstIndex + 4);
2271 
2272  Mat3x3 RDelta(CGR_Rot::MatR, gCurr);
2273 
2274  RCurr = RDelta*RRef;
2275  WCurr = Vec3(X, iFirstIndex + 10);
2276 }
2277 
2278 /* Inverse Dynamics: */
2279 void
2280 StructNode::Update(const VectorHandler& X, InverseDynamics::Order iOrder)
2281 {
2282 #ifdef USE_AUTODIFF
2283  bUpdateRotation = true;
2284 #endif
2285 
2286  integer iFirstIndex = iGetFirstIndex();
2287  switch (iOrder) {
2288  case InverseDynamics::POSITION: {
2289  XCurr = Vec3(X, iFirstIndex + 1);
2290  gCurr = Vec3(X, iFirstIndex + 4);
2291  Mat3x3 RDelta(CGR_Rot::MatR, gCurr);
2292  RCurr = RDelta*RRef;
2293  // gCurr = ::Zero3;
2294  } break;
2295 
2296  case InverseDynamics::VELOCITY: {
2297  VCurr = Vec3(X, iFirstIndex + 1);
2298 #if 0
2299  gPCurr = Vec3(X, iFirstIndex + 4);
2300  Mat3x3 RDelta(CGR_Rot::MatR, gCurr);
2301  WCurr = Mat3x3(CGR_Rot::MatG, gCurr)*gPCurr + RDelta*WRef;
2302 #endif
2303  WCurr = Vec3(X, iFirstIndex + 4);
2304  } break;
2305 
2306  case InverseDynamics::ACCELERATION: {
2307  XPPCurr = Vec3(X, iFirstIndex + 1);
2308  WPCurr = Vec3(X, iFirstIndex + 4);
2309  } break;
2310 
2311  default:
2312  ASSERT(0);
2313  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2314  }
2315 }
2316 
2317 /* Funzioni di inizializzazione, ereditate da DofOwnerOwner */
2318 void
2319 StructNode::SetInitialValue(VectorHandler& X)
2320 {
2321  /* FIXME: why is this called? */
2322  integer iIndex = iGetFirstIndex();
2323 
2324  X.Put(iIndex + 1, XCurr);
2325  X.Put(iIndex + 4, Zero3);
2326  X.Put(iIndex + 7, VCurr);
2327  X.Put(iIndex + 10, WCurr);
2328 }
2329 
2330 
2331 void
2332 StructNode::SetValue(DataManager *pDM,
2333  VectorHandler& X, VectorHandler& XP,
2334  SimulationEntity::Hints *ph)
2335 {
2336 #ifdef USE_AUTODIFF
2337  bUpdateRotation = true;
2338 #endif
2339 
2340 #ifdef MBDYN_X_RELATIVE_PREDICTION
2341  if (pRefNode) {
2342  Vec3 Xtmp = XPrev - pRefNode->GetXCurr();
2343  const Mat3x3& R0 = pRefNode->GetRCurr();
2344  const Vec3& V0 = pRefNode->GetVCurr();
2345  const Vec3& W0 = pRefNode->GetWCurr();
2346 
2347  XPrev = R0.MulTV(Xtmp);
2348  RPrev = R0.MulTM(RCurr);
2349  VPrev = R0.MulTV(VCurr - V0 - W0.Cross(Xtmp));
2350  WPrev = R0.MulTV(WCurr - W0);
2351 
2352 #if 0
2353  std::cout << "StructNode(" << GetLabel() << "): "
2354  "SetValue: X=" << XPrev
2355  << ", R=" << RPrev
2356  << ", V=" << VPrev
2357  << ", W=" << WPrev
2358  << std::endl;
2359 #endif
2360 
2361  } else
2362 #endif /* MBDYN_X_RELATIVE_PREDICTION */
2363  {
2364  /* FIXME: in any case, we start with Crank-Nicolson ... */
2365  XPrev = XCurr;
2366  RPrev = RCurr;
2367  VPrev = VCurr;
2368  WPrev = WCurr;
2369  // Without the next line the Jacobian matrix of most elements
2370  // will be incorrect during the initial derivatives calculation
2371  // if RCurr has been changed during initial assembly!
2372  // The reason is, that most elements assume that
2373  // RCurr = RDelta * RRef
2374  // and not RCurr = RDelta * RPrev
2375  RRef = RCurr;
2376  WRef = WCurr;
2377  }
2378 
2379  integer iFirstIndex = iGetFirstIndex();
2380  X.Put(iFirstIndex + 1, XPrev);
2381  X.Put(iFirstIndex + 4, Zero3);
2382  gRef = gCurr = gPRef = gPCurr = Zero3;
2383  XP.Put(iFirstIndex + 1, VPrev);
2384  XP.Put(iFirstIndex + 4, WPrev);
2385 }
2386 
2387 
2388 void
2389 StructNode::BeforePredict(VectorHandler& X,
2390  VectorHandler& XP,
2391  VectorHandler& XPr,
2392  VectorHandler& XPPr) const
2393 {
2394  integer iFirstPos = iGetFirstIndex();
2395 
2396 #ifdef MBDYN_X_RELATIVE_PREDICTION
2397  /* If pRefNode is defined, the prediction is made
2398  * on the data in the reference frame it provides */
2399  if (pRefNode) {
2400 
2401  /*
2402  x_r = R_0^T * ( x - x_0 )
2403  R_r = R_0^T * R
2404  v_r = R_0^T * ( v - v_0 - omega_0 \times ( x - x_0 ) )
2405  omega_r = R_0^T * ( omega - omega_0 )
2406  */
2407  Vec3 Xtmp = XCurr - pRefNode->GetXCurr();
2408  const Mat3x3& R0 = pRefNode->GetRCurr();
2409  const Vec3& V0 = pRefNode->GetVCurr();
2410  const Vec3& W0 = pRefNode->GetWCurr();
2411 
2412  XCurr = R0.MulTV(Xtmp);
2413  RCurr = R0.MulTM(RCurr);
2414  VCurr = R0.MulTV(VCurr - V0 - W0.Cross(Xtmp));
2415  WCurr = R0.MulTV(WCurr - W0);
2416 
2417  /* update state vectors with relative position and velocity */
2418  X.Put(iFirstPos + 1, XCurr);
2419  XP.Put(iFirstPos + 1, VCurr);
2420  XPr.Put(iFirstPos + 1, XPrev);
2421  XPPr.Put(iFirstPos + 1, VPrev);
2422 
2423 #if 0
2424  std::cout << "StructNode(" << GetLabel() << "): "
2425  "BeforePredict: X=" << XCurr
2426  << ", R=" << RCurr
2427  << ", V=" << VCurr
2428  << ", W=" << WCurr
2429  << std::endl;
2430 #endif
2431  }
2432 #endif /* MBDYN_X_RELATIVE_PREDICTION */
2433 
2434  /* Questa e' la predizione "consistente", ovvero usa come gdl
2435  * di rotazione i parametri di rotazione "totali" per predire
2436  * la configurazione al nuovo passo, quindi ritorna in forma
2437  * incrementale */
2438 
2439  /* Calcolo la matrice RDelta riferita a tutto il passo trascorso
2440  * all'indietro */
2441  Mat3x3 RDelta(RPrev.MulMT(RCurr));
2442 
2443  /* Mi assicuro che g al passo corrente sia nullo */
2444  X.Put(iFirstPos + 4, Zero3);
2445 
2446  /* Calcolo g al passo precedente attraverso la matrice RDelta riferita
2447  * a tutto il passo. Siccome RDelta e' calcolata all'indietro,
2448  * i parametri sono gia' con il segno corretto */
2449  Vec3 gPrev(CGR_Rot::Param, RDelta);
2450  XPr.Put(iFirstPos + 4, gPrev);
2451 
2452  /* Calcolo gP al passo precedente attraverso la definizione
2453  * mediante le Omega. Siccome i parametri sono con il segno meno
2454  * e la matrice RDelta e' gia' calcolata all'indietro, l'insieme
2455  * e' consistente */
2456  XPPr.Put(iFirstPos + 4, Mat3x3(CGR_Rot::MatGm1, gPrev)*WPrev);
2457 
2458  /* Metto Omega al passo corrente come gP (perche' G(0) = I) */
2459  XP.Put(iFirstPos + 4, WCurr);
2460 
2461 #if 0
2462  std::cout
2463  << " " << std::setw(16) << "prev" << std::setw(16) << "curr" << std::setw(16) << GetLabel() << std::endl
2464  << "x:" << std::setw(16) << XPrev(1) << std::setw(16) << XCurr(1) << std::endl
2465  << " " << std::setw(16) << XPrev(2) << std::setw(16) << XCurr(2) << std::endl
2466  << " " << std::setw(16) << XPrev(3) << std::setw(16) << XCurr(3) << std::endl
2467  << "v:" << std::setw(16) << VPrev(1) << std::setw(16) << VCurr(1) << std::endl
2468  << " " << std::setw(16) << VPrev(2) << std::setw(16) << VCurr(2) << std::endl
2469  << " " << std::setw(16) << VPrev(3) << std::setw(16) << VCurr(3) << std::endl
2470  << "g:" << std::setw(16) << gPrev(1) << std::setw(16) << 0 << std::endl
2471  << " " << std::setw(16) << gPrev(2) << std::setw(16) << 0 << std::endl
2472  << " " << std::setw(16) << gPrev(3) << std::setw(16) << 0 << std::endl
2473  << "w:" << std::setw(16) << XP(iFirstPos+4) << std::setw(16) << WCurr(1) << std::endl
2474  << " " << std::setw(16) << XP(iFirstPos+5) << std::setw(16) << WCurr(2) << std::endl
2475  << " " << std::setw(16) << XP(iFirstPos+6) << std::setw(16) << WCurr(3) << std::endl;
2476 #endif
2477 
2478  XPrev = XCurr;
2479  VPrev = VCurr;
2480 
2481  /* Pongo la R al passo precedente uguale a quella corrente
2482  * mi servira' se devo ripetere il passo con un diverso Delta t
2483  * e per la rettifica dopo la predizione */
2484  RPrev = RCurr;
2485 
2486  /* Pongo le Omega al passo precedente uguali alle Omega al passo corrente
2487  * mi servira' per la correzione dopo la predizione */
2488  WPrev = WCurr;
2489 }
2490 
2491 void
2492 StructNode::AfterPredict(VectorHandler& X, VectorHandler& XP)
2493 {
2494 #ifdef USE_AUTODIFF
2495  bUpdateRotation = true;
2496 #endif
2497 
2498  integer iFirstIndex = iGetFirstIndex();
2499 
2500  /* Spostamento e velocita' aggiornati */
2501  XCurr = Vec3(X, iFirstIndex + 1);
2502  VCurr = Vec3(XP, iFirstIndex + 1);
2503 
2504  /* Ottengo il g predetto */
2505  gRef = Vec3(X, iFirstIndex + 4);
2506 
2507  /* Calcolo la matrice RDelta derivante dalla predizione */
2508  Mat3x3 RDelta(CGR_Rot::MatR, gRef);
2509 
2510  /* Calcolo la R corrente in base alla predizione */
2511  RCurr = RDelta*RPrev;
2512 
2513  /* Calcolo la Omega corrente in base alla predizione (gP "totale") */
2514  gPRef = Vec3(XP, iFirstIndex + 4);
2515 
2516  /* Calcolo il nuovo Omega */
2517  WCurr = Mat3x3(CGR_Rot::MatG, gRef)*gPRef;
2518 
2519  /* Resetto i parametri di rotazione e le derivate, g e gP */
2520  X.Put(iFirstIndex + 4, Zero3);
2521  XP.Put(iFirstIndex + 4, Zero3);
2522 
2523  gCurr = gPCurr = Zero3;
2524 
2525 #ifdef MBDYN_X_RELATIVE_PREDICTION
2526  if (pRefNode) {
2527 
2528  /*
2529  x = x_0 + R_0 * x_r
2530  R = R_0 * R_r
2531  v = v_0 + omega_0 \times ( R_0 * x_r ) + R_0 * v_r
2532  omega = omega_0 + R_0 * omega_r
2533  */
2534  Vec3 X0 = pRefNode->GetXCurr();
2535  Mat3x3 R0 = pRefNode->GetRCurr();
2536  Vec3 V0 = pRefNode->GetVCurr();
2537  Vec3 W0 = pRefNode->GetWCurr();
2538 
2539  XCurr = R0*XCurr; /* temporary */
2540  RCurr = R0*RCurr;
2541  VCurr = V0 + W0.Cross(XCurr) + R0*VCurr;
2542  WCurr = W0 + R0*WCurr;
2543  XCurr += X0; /* plus reference */
2544 
2545  /* alcuni usano anche le predizioni dei parametri
2546  * di rotazione e delle loro derivate come riferimento
2547  * (approccio updated-updated); quindi calcolo
2548  * i parametri di riferimento come i parametri
2549  * che danno una predizione pari alla variazione
2550  * di R0 piu' l'incremento relativo, e le derivate
2551  * dei parametri corrispondenti */
2552  gRef = Vec3(CGR_Rot::Param, R0*RDelta.MulMT(pRefNode->GetRPrev()));
2553  gPRef = Mat3x3(CGR_Rot::MatGm1, gRef)*WCurr;
2554 
2555  /* to be safe, the correct values are put back
2556  * in the state vectors */
2557  X.Put(iFirstIndex + 1, XCurr);
2558  XP.Put(iFirstIndex + 1, VCurr);
2559 
2560 #if 0
2561  std::cout << "StructNode(" << GetLabel() << "): "
2562  "AfterPredict: X=" << XCurr
2563  << ", R=" << RCurr
2564  << ", V=" << VCurr
2565  << ", W=" << WCurr
2566  << std::endl;
2567 #endif
2568  }
2569 #endif /* MBDYN_X_RELATIVE_PREDICTION */
2570 
2571  RRef = RCurr;
2572  WRef = WCurr;
2573 
2574 #if 0
2575  /* Ortho check */
2576  Mat3x3 RRT = RCurr.MulTM(RCurr);
2577  RRT(1, 1) -= 1.;
2578  RRT(2, 2) -= 1.;
2579  RRT(3, 3) -= 1.;
2580  doublereal dmax = 0.;
2581  for (int r = 1; r <= 3; r++) {
2582  for (int c = 1; c <= 3; c++) {
2583  dmax = std::max(dmax, fabs(RRT(r, c)));
2584  }
2585  }
2586  silent_cout("### StructNode(" << GetLabel() << ") " << dmax << std::endl);
2587 #endif
2588 }
2589 
2590 /* Inverse Dynamics: */
2591 void
2592 StructNode::AfterConvergence(const VectorHandler& X,
2593  const VectorHandler& XP,
2594  const VectorHandler& XPP)
2595 {
2596 #ifdef USE_AUTODIFF
2597  bUpdateRotation = true;
2598 #endif
2599 /* Right now, AfterConvergence is performed only on position
2600  * to reset orientation parameters. XPrime and XPrimePrime are
2601  * left for compatibility with the virtual method in
2602  * class SimulationEntity */
2603 
2604  integer iFirstIndex = iGetFirstIndex();
2605 
2606 
2607  /* Orientation Parameters:
2608  * Get g and impose it as gRef: successive iterations
2609  * use gRef as reference and the solution is a perturbation
2610  * from it */
2611  gRef = Vec3(X, iFirstIndex + 4);
2612  gCurr = Zero3;
2613  RRef = RCurr;
2614  WRef = WCurr;
2615 
2616  XPrev = XCurr;
2617  RPrev = RCurr;
2618  VPrev = VCurr;
2619  WPrev = WCurr;
2620  XPPPrev = XPPCurr;
2621  WPPrev = WPCurr;
2622 }
2623 
2624 /*
2625  * Metodi per l'estrazione di dati "privati".
2626  * Si suppone che l'estrattore li sappia interpretare.
2627  * Come default non ci sono dati privati estraibili
2628  */
2629 unsigned int
2630 StructNode::iGetNumPrivData(void) const
2631 {
2632  unsigned i =
2633  3 // X
2634  + 3 // x (R^T * X)
2635  + 3 // Phi
2636  + 3 // XP
2637  + 3 // xP (R^T * XP)
2638  + 3 // Omega
2639  + 3 // omega (R^T * Omega)
2640  + 3 // Euler angles (123)
2641  + 3 // Euler angles (313)
2642  + 3 // Euler angles (321)
2643  + 4; // Euler parameters
2644 
2645  if (bComputeAccelerations()) {
2646  i +=
2647  3 // XPP
2648  + 3 // xPP (R^T * XPP)
2649  + 3 // OmegaP
2650  + 3; // omegaP (R^T * OmegaP)
2651  }
2652 
2653  return i;
2654 }
2655 
2656 /*
2657  * Maps a string (possibly with substrings) to a private data;
2658  * returns a valid index ( > 0 && <= iGetNumPrivData()) or 0
2659  * in case of unrecognized data; error must be handled by caller
2660  */
2661 unsigned int
2662 StructNode::iGetPrivDataIdx(const char *s) const
2663 {
2664  long idx;
2665  char *next;
2666  std::string sDataName(s);
2667 
2668  const char *brk = std::strchr(s, '[' /*]*/ );
2669  if (brk == 0) {
2670  return 0;
2671  }
2672 
2673  size_t len = brk - s;;
2674  brk++;
2675 
2676  errno = 0;
2677  idx = strtol(brk, &next, 10);
2678  int save_errno = errno;
2679  if (next == brk || strcmp(next, /*[*/ "]") != 0) {
2680  return 0;
2681  }
2682 
2683  if (save_errno == ERANGE) {
2684  silent_cerr("StructNode(" << GetLabel() << "): "
2685  "warning, private data index "
2686  << std::string(brk, next - brk)
2687  << " overflows" << std::endl);
2688  return 0;
2689  }
2690 
2691  /*
2692  X 0 + idx idx = {1,3}
2693  x 3 + idx idx = {1,3}
2694  Phi 6 + idx idx = {1,3}
2695  XP 9 + idx idx = {1,3}
2696  x 12 + idx idx = {1,3}
2697  Omega 15 + idx idx = {1,3}
2698  omega 18 + idx idx = {1,3}
2699  E | E123 21 + idx idx = {1,3}
2700  E313 24 + idx idx = {1,3}
2701  E321 27 + idx idx = {1,3}
2702  PE 31 + idx idx = {0,3}
2703  -------------------------------------------
2704  XPP 34 + idx idx = {1,3}
2705  xPP 37 + idx idx = {1,3}
2706  OmegaP 40 + idx idx = {1,3}
2707  omegaP 43 + idx idx = {1,3}
2708  */
2709 
2710  if (strncmp(s, "PE", len) == 0) {
2711  if (idx < 0 || idx > 3) {
2712  return 0;
2713  }
2714 
2715  return 31 + idx;
2716  }
2717 
2718  if (idx < 1 || idx > 3) {
2719  return 0;
2720  }
2721 
2722  if (strncmp(s, "X", len) == 0) {
2723  return 0 + idx;
2724  }
2725 
2726  if (strncmp(s, "x", len) == 0) {
2727  return 3 + idx;
2728  }
2729 
2730  if (strncmp(s, "Phi", len) == 0) {
2731  return 6 + idx;
2732  }
2733 
2734  if (strncmp(s, "XP", len) == 0) {
2735  return 9 + idx;
2736  }
2737 
2738  if (strncmp(s, "xP", len) == 0) {
2739  return 12 + idx;
2740  }
2741 
2742  if (strncmp(s, "Omega", len) == 0) {
2743  return 15 + idx;
2744  }
2745 
2746  if (strncmp(s, "omega", len) == 0) {
2747  return 18 + idx;
2748  }
2749 
2750  if (strncmp(s, "E", len) == 0
2751  || strncmp(s, "E123", len) == 0)
2752  {
2753  return 21 + idx;
2754  }
2755 
2756  if (strncmp(s, "E313", len) == 0) {
2757  return 24 + idx;
2758  }
2759 
2760  if (strncmp(s, "E321", len) == 0) {
2761  return 27 + idx;
2762  }
2763 
2764  bool bca = false;
2765  unsigned i;
2766  if (strncmp(s, "XPP", len) == 0) {
2767  bca = true;
2768  i = 34 + idx;
2769 
2770  } else if (strncmp(s, "xPP", len) == 0) {
2771  bca = true;
2772  i = 37 + idx;
2773 
2774  } else if (strncmp(s, "OmegaP", len) == 0) {
2775  bca = true;
2776  i = 40 + idx;
2777 
2778  } else if (strncmp(s, "omegaP", len) == 0) {
2779  bca = true;
2780  i = 43 + idx;
2781 
2782  } else {
2783  // error
2784  return 0;
2785  }
2786 
2787  // NOTE: bComputeAccels is set only if iGetPrivDataIdx() is called
2788  // first; it is not when the (deprecated) idx is directly used.
2789  if (bca) {
2790  if (!const_cast<StructNode *>(this)->ComputeAccelerations(true)) {
2791  silent_cerr("StructNode(" << GetLabel() << "): "
2792  "request to compute accelerations failed, requested by private data \"" << sDataName << "\""
2793  << std::endl);
2794  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2795  }
2796  }
2797 
2798  return i;
2799 }
2800 
2801 /*
2802  * Returns the current value of a private data
2803  * with 0 < i <= iGetNumPrivData()
2804  */
2805 doublereal
2806 StructNode::dGetPrivData(unsigned int i) const
2807 {
2808  switch (i) {
2809  case 1:
2810  case 2:
2811  case 3:
2812  return XCurr(i);
2813 
2814  case 4:
2815  case 5:
2816  case 6:
2817  return RCurr.GetVec(i - 3)*XCurr;
2818 
2819  case 7:
2820  case 8:
2821  case 9: {
2822  /* TODO */
2823  Vec3 Phi(RotManip::VecRot(RCurr));
2824  return Phi(i - 6);
2825  }
2826 
2827  case 10:
2828  case 11:
2829  case 12:
2830  return VCurr(i - 9);
2831 
2832  case 13:
2833  case 14:
2834  case 15:
2835  return RCurr.GetVec(i - 12)*VCurr;
2836 
2837  case 16:
2838  case 17:
2839  case 18:
2840  return WCurr(i - 15);
2841 
2842  case 19:
2843  case 20:
2844  case 21:
2845  return RCurr.GetVec(i - 18)*WCurr;
2846 
2847  case 22:
2848  case 23:
2849  case 24: {
2850  Vec3 Phi(MatR2EulerAngles123(RCurr));
2851  return Phi(i - 21);
2852  }
2853 
2854  case 25:
2855  case 26:
2856  case 27: {
2857  Vec3 Phi(MatR2EulerAngles313(RCurr));
2858  return Phi(i - 24);
2859  }
2860 
2861  case 28:
2862  case 29:
2863  case 30: {
2864  Vec3 Phi(MatR2EulerAngles321(RCurr));
2865  return Phi(i - 27);
2866  }
2867 
2868  case 31:
2869  case 32:
2870  case 33:
2871  case 34: {
2872  /* TODO */
2873  Vec3 e;
2874  doublereal e0;
2875  MatR2EulerParams(RCurr, e0, e);
2876  if (i == 31) {
2877  return e0;
2878  }
2879  return e(i - 31);
2880  }
2881 
2882  case 35:
2883  case 36:
2884  case 37:
2885  ASSERT(bComputeAccelerations() == true);
2886  return XPPCurr(i - 34);
2887 
2888  case 38:
2889  case 39:
2890  case 40:
2891  ASSERT(bComputeAccelerations() == true);
2892  return RCurr.GetVec(i - 37)*XPPCurr;
2893 
2894  case 41:
2895  case 42:
2896  case 43:
2897  ASSERT(bComputeAccelerations() == true);
2898  return WPCurr(i - 40);
2899 
2900  case 44:
2901  case 45:
2902  case 46:
2903  ASSERT(bComputeAccelerations() == true);
2904  return RCurr.GetVec(i - 43)*WPCurr;
2905  }
2906 
2907  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2908 }
2909 
2910 /* StructNode - end */
2911 
2912 
2913 /* DynamicStructNode - begin */
2914 
2915 DynamicStructNode::DynamicStructNode(unsigned int uL,
2916  const DofOwner* pDO,
2917  const Vec3& X0,
2918  const Mat3x3& R0,
2919  const Vec3& V0,
2920  const Vec3& W0,
2921  const StructNode *pRN,
2922  const RigidBodyKinematics *pRBK,
2923  doublereal dPosStiff,
2924  doublereal dVelStiff,
2925  bool bOmRot,
2926  OrientationDescription ood,
2927  flag fOut)
2928 :
2929 StructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut),
2930 DynamicStructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut),
2931 StructNode(uL, pDO, X0, R0, V0, W0, pRN, pRBK, dPosStiff, dVelStiff, bOmRot, ood, fOut)
2932 {
2933  NO_OP;
2934 }
2935 
2936 
2937 /* Distruttore (per ora e' banale) */
2938 DynamicStructNode::~DynamicStructNode(void)
2939 {
2940  NO_OP;
2941 }
2942 
2943 
2944 /* Tipo di nodo strutturale */
2945 StructNode::Type
2946 DynamicStructNode::GetStructNodeType(void) const
2947 {
2948  return StructNode::DYNAMIC;
2949 }
2950 
2951 const Vec3&
2952 DynamicStructNode::GetWP(void) const
2953 {
2954  return GetWPCurr();
2955 }
2956 
2957 std::ostream&
2958 DynamicStructNode::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
2959 {
2960  integer iIndex = iGetFirstIndex();
2961 
2962  StructNode::DescribeDof(out, prefix, bInitial);
2963 
2964  if (bInitial == false) {
2965  out
2966  << prefix << iIndex + 7 << "->" << iIndex + 9 << ": "
2967  "momentum [Bx,By,Bz]" << std::endl
2968  << prefix << iIndex + 10 << "->" << iIndex + 12 << ": "
2969  "momenta moment [Gx,Gy,Gz]" << std::endl;
2970  }
2971 
2972  return out;
2973 }
2974 
2975 void
2976 DynamicStructNode::DescribeDof(std::vector<std::string>& desc, bool bInitial, int i) const
2977 {
2978  if (bInitial || i == -1 || (i >= 0 && i < 6)) {
2979  StructNode::DescribeDof(desc, bInitial, i);
2980 
2981  if (bInitial || (i >= 0 && i < 6)) {
2982  return;
2983  }
2984  }
2985 
2986  if (i == -1) {
2987  desc.resize(12);
2988 
2989  } else {
2990  desc.resize(1);
2991  }
2992 
2993  std::ostringstream os;
2994  os << "StructNode(" << GetLabel() << ")";
2995 
2996  if (i == -1) {
2997  std::string name = os.str();
2998 
2999  for (i = 6; i < 12; i++) {
3000  os.str(name);
3001  os.seekp(0, std::ios_base::end);
3002  os << ": " << sn_dof[i/3] << xyz[i%3];
3003  desc[i] = os.str();
3004  }
3005 
3006  } else {
3007  if (i < 6 || i >= 12) {
3008  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3009  }
3010 
3011  os << ": " << sn_dof[i/3] << xyz[i%3];
3012  desc[0] = os.str();
3013  }
3014 }
3015 
3016 std::ostream&
3017 DynamicStructNode::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
3018 {
3019  if (bInitial == false) {
3020  integer iIndex = iGetFirstIndex();
3021 
3022  out
3023  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
3024  "momentum definition [Bx,By,Bz]" << std::endl
3025  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
3026  "momenta moment definition [Gx,Gy,Gz]" << std::endl;
3027  }
3028 
3029  StructNode::DescribeEq(out, prefix, bInitial);
3030 
3031  return out;
3032 }
3033 
3034 void
3035 DynamicStructNode::DescribeEq(std::vector<std::string>& desc, bool bInitial, int i) const
3036 {
3037  if (bInitial || i == -1 || (i >= 6 && i < 12)) {
3038  int new_i = i;
3039  if (!bInitial && i != -1) {
3040  new_i = i - 6;
3041  }
3042  StructNode::DescribeEq(desc, bInitial, new_i);
3043 
3044  if (bInitial || (i >= 6 && i < 12)) {
3045  return;
3046  }
3047  }
3048 
3049  if (i == -1) {
3050  desc.resize(12);
3051  for (int j = 0; j < 6; j++) {
3052  desc[6 + j] = desc[j];
3053  }
3054 
3055  } else {
3056  desc.resize(1);
3057  }
3058 
3059  std::ostringstream os;
3060  os << "StructNode(" << GetLabel() << ")";
3061 
3062  if (i == -1) {
3063  std::string name(os.str());
3064 
3065  for (i = 0; i < 6; i++) {
3066  os.str(name);
3067  os.seekp(0, std::ios_base::end);
3068  os << ": " << sn_eq[i/3] << xyz[i%3];
3069  desc[i] = os.str();
3070  }
3071 
3072  } else {
3073  if (i < 0 || i >= 6) {
3074  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3075  }
3076 
3077  os << ": " << sn_eq[i/3] << xyz[i%3];
3078  desc[0] = os.str();
3079  }
3080 }
3081 
3082 /* Usato dalle forze astratte, dai bulk ecc., per assemblare le forze
3083  * al posto giusto */
3084 integer
3085 DynamicStructNode::iGetFirstRowIndex(void) const
3086 {
3087  return iGetFirstMomentumIndex();
3088 }
3089 
3090 /* delegate to autostr node */
3091 void
3092 DynamicStructNode::AddInertia(const doublereal& dm, const Vec3& dS,
3093  const Mat3x3& dJ) const
3094 {
3095  /* FIXME: do it only if to be output... */
3096  if (bComputeAccelerations()) {
3097  dynamic_cast<AutomaticStructElem *>(pAutoStr)->AddInertia(dm, dS, dJ);
3098  }
3099 }
3100 
3101 /* Accesso ai suoi dati */
3102 const Vec3&
3103 DynamicStructNode::GetGCurr(void) const
3104 {
3105  return pAutoStr->GetGCurr();
3106 }
3107 
3108 const Vec3&
3109 DynamicStructNode::GetGPCurr(void) const
3110 {
3111  return pAutoStr->GetGPCurr();
3112 }
3113 
3114 void
3115 DynamicStructNode::Update(const VectorHandler& X, const VectorHandler& XP)
3116 {
3117  StructNode::Update(X, XP);
3118  if (bComputeAccelerations()) {
3119  /* FIXME: pAutoStr is 0 in ModalNode */
3120  ASSERT(pAutoStr != 0);
3121 
3122  // FIXME: based on values set during previous
3123  // of AutomaticStructural::AssRes()
3124  dynamic_cast<const AutomaticStructElem *>(pAutoStr)->ComputeAccelerations(XPPCurr, WPCurr);
3125  }
3126 }
3127 
3128 void
3129 DynamicStructNode::AfterConvergence(const VectorHandler& X,
3130  const VectorHandler& XP)
3131 {
3132  if (bComputeAccelerations()) {
3133  /* FIXME: pAutoStr is 0 in ModalNode */
3134  ASSERT(pAutoStr != 0);
3135 
3136  // FIXME: based on values set during previous
3137  // of AutomaticStructural::AssRes()
3138  dynamic_cast<const AutomaticStructElem *>(pAutoStr)->ComputeAccelerations(XPPCurr, WPCurr);
3139  }
3140 }
3141 
3142 void
3143 DynamicStructNode::BeforePredict(VectorHandler& X,
3144  VectorHandler& XP,
3145  VectorHandler& XPr,
3146  VectorHandler& XPPr) const
3147 {
3148  if (bComputeAccelerations()) {
3149  XPPPrev = XPPCurr;
3150  WPPrev = WPCurr;
3151  }
3152 
3153  StructNode::BeforePredict(X, XP, XPr, XPPr);
3154 }
3155 
3156 /* Restituisce il valore del dof iDof;
3157  * se differenziale, iOrder puo' essere = 1 per la derivata */
3158 const doublereal&
3159 DynamicStructNode::dGetDofValue(int iDof, int iOrder) const
3160 {
3161  ASSERT(iDof >= 1 && iDof <= 6);
3162  ASSERT(iOrder >= 0 && iOrder <= 2);
3163 
3164  if (iOrder == 2) {
3165  /* FIXME: should not happen */
3166  ASSERT(bComputeAccelerations());
3167  if (!bComputeAccelerations()) {
3168  silent_cerr("DynamicStructNode::dGetDofValue("
3169  << iDof << "," << iOrder << "): "
3170  "accelerations are not computed while they should"
3171  << std::endl);
3172  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3173  }
3174 
3175 #if 1
3176  /* FIXME: might need to compute them in order to be
3177  * as up to date as possible; however, elements that contribute
3178  * to inertia should assemble first...
3179  */
3180  dynamic_cast<const AutomaticStructElem *>(pAutoStr)->ComputeAccelerations(XPPCurr, WPCurr);
3181 #endif
3182 
3183  if (iDof >= 1 && iDof <= 3) {
3184  return XPPCurr(iDof);
3185  } else {
3186  return WPCurr(iDof - 3);
3187  }
3188  }
3189 
3190  return StructNode::dGetDofValue(iDof, iOrder);
3191 }
3192 
3193 /* Restituisce il valore del dof iDof al passo precedente;
3194  * se differenziale, iOrder puo' essere = 1 per la derivata */
3195 const doublereal&
3196 DynamicStructNode::dGetDofValuePrev(int iDof, int iOrder) const
3197 {
3198  ASSERT(iDof >= 1 && iDof <= 6);
3199  ASSERT(iOrder == 0 || iOrder == 1);
3200 
3201  if (iOrder == 2) {
3202  /* FIXME: should not happen */
3203  ASSERT(bComputeAccelerations());
3204  if (!bComputeAccelerations()) {
3205  silent_cerr("DynamicStructNode::dGetDofValuePrev("
3206  << iDof << "," << iOrder << "): "
3207  "accelerations are not computed while they should"
3208  << std::endl);
3209  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3210  }
3211 
3212  if (iDof >= 1 && iDof <= 3) {
3213  return XPPPrev(iDof);
3214  } else {
3215  return WPPrev(iDof - 3);
3216  }
3217  }
3218 
3219  return StructNode::dGetDofValuePrev(iDof, iOrder);
3220 }
3221 
3222 /* Setta il valore del dof iDof a dValue;
3223  * se differenziale, iOrder puo' essere = 1 per la derivata */
3224 void
3225 DynamicStructNode::SetDofValue(const doublereal& dValue,
3226  unsigned int iDof,
3227  unsigned int iOrder /* = 0 */ )
3228 {
3229  ASSERT(iDof >= 1 && iDof <= 6);
3230  ASSERT(iOrder == 0 || iOrder == 1);
3231 
3232  if (iOrder == 2) {
3233  /* FIXME: should not happen */
3234  ASSERT(bComputeAccelerations());
3235  if (!bComputeAccelerations()) {
3236  silent_cerr("DynamicStructNode::SetDofValue("
3237  << dValue << "," << iDof << "," << iOrder << "): "
3238  "accelerations are not computed while they should"
3239  << std::endl);
3240  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3241  }
3242 
3243  if (iDof >= 1 && iDof <= 3) {
3244  XPPCurr(iDof) = dValue;
3245 
3246  } else {
3247  WPCurr(iDof - 3) = dValue;
3248  }
3249 
3250  } else {
3251  StructNode::SetDofValue(iDof, iOrder);
3252  }
3253 }
3254 
3255 /* DynamicStructNode - end */
3256 
3257 
3258 /* StaticStructNode - begin */
3259 
3260 /* Costruttore definitivo */
3261 StaticStructNode::StaticStructNode(unsigned int uL,
3262  const DofOwner* pDO,
3263  const Vec3& X0,
3264  const Mat3x3& R0,
3265  const Vec3& V0,
3266  const Vec3& W0,
3267  const StructNode *pRN,
3268  const RigidBodyKinematics *pRBK,
3269  doublereal dPosStiff,
3270  doublereal dVelStiff,
3271  bool bOmRot,
3272  OrientationDescription ood,
3273  flag fOut)
3274 :
3275 StructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut),
3276 StaticStructDispNode(uL, pDO, X0, V0, pRN, pRBK, dPosStiff, dVelStiff, ood, fOut),
3277 StructNode(uL, pDO, X0, R0, V0, W0, pRN, pRBK, dPosStiff, dVelStiff, bOmRot,
3278  ood, fOut)
3279 {
3280  NO_OP;
3281 }
3282 
3283 
3284 /* Distruttore (per ora e' banale) */
3285 StaticStructNode::~StaticStructNode(void)
3286 {
3287  NO_OP;
3288 }
3289 
3290 
3291 /* Tipo di nodo strutturale */
3292 StructNode::Type
3293 StaticStructNode::GetStructNodeType(void) const
3294 {
3295  return StructNode::STATIC;
3296 }
3297 
3298 /* StaticStructNode - end */
3299 
3300 
3301 /* ModalNode - begin */
3302 
3303 ModalNode::ModalNode(unsigned int uL,
3304  const DofOwner* pDO,
3305  const Vec3& X0,
3306  const Mat3x3& R0,
3307  const Vec3& V0,
3308  const Vec3& W0,
3309  const RigidBodyKinematics *pRBK,
3310  doublereal dPosStiff,
3311  doublereal dVelStiff,
3312  bool bOmRot,
3313  OrientationDescription ood,
3314  flag fOut)
3315 :
3316 StructDispNode(uL, pDO, X0, V0, 0, pRBK, dPosStiff, dVelStiff, ood, fOut),
3317 DynamicStructNode(uL, pDO, X0, R0, V0, W0, 0, pRBK,
3318  dPosStiff, dVelStiff, bOmRot, ood, fOut)
3319 {
3320  /* XPP and WP are not known in ModalNode */
3321  ComputeAccelerations(false);
3322 }
3323 
3324 
3325 /* Distruttore (per ora e' banale) */
3326 ModalNode::~ModalNode(void)
3327 {
3328  NO_OP;
3329 }
3330 
3331 
3332 /* Tipo di nodo strutturale */
3333 StructNode::Type
3334 ModalNode::GetStructNodeType(void) const
3335 {
3336  return StructNode::MODAL;
3337 }
3338 
3339 
3340 /* Usato dalle forze astratte, dai bulk ecc., per assemblare le forze
3341  * al posto giusto */
3342 integer
3343 ModalNode::iGetFirstRowIndex(void) const
3344 {
3345  return iGetFirstMomentumIndex();
3346 }
3347 
3348 std::ostream&
3349 ModalNode::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
3350 {
3351  StructNode::DescribeDof(out, prefix, bInitial);
3352 
3353  if (bInitial == false) {
3354  integer iIndex = iGetFirstIndex();
3355 
3356  out
3357  << prefix << iIndex + 7 << "->" << iIndex + 9 << ": "
3358  "velocity [vx,vy,vz]" << std::endl
3359  << prefix << iIndex + 10 << "->" << iIndex + 12 << ": "
3360  "angular velocity [wx,wy,wz]" << std::endl;
3361  }
3362 
3363  return out;
3364 }
3365 
3366 void
3367 ModalNode::DescribeDof(std::vector<std::string>& desc, bool bInitial, int i) const
3368 {
3369  if (bInitial || i == -1 || (i >= 0 && i < 6)) {
3370  StructNode::DescribeDof(desc, bInitial, i);
3371 
3372  if (bInitial || (i >= 0 && i < 6)) {
3373  for (size_t ii = 0; ii < desc.size(); ii++) {
3374  desc[ii] = "Modal" + desc[ii];
3375  }
3376  return;
3377  }
3378  }
3379 
3380  if (i == -1) {
3381  desc.resize(12);
3382 
3383  } else {
3384  desc.resize(1);
3385  }
3386 
3387  std::ostringstream os;
3388  os << "ModalStructNode(" << GetLabel() << ")";
3389 
3390  if (i == -1) {
3391  std::string name = os.str();
3392 
3393  for (i = 0; i < 6; i++) {
3394  desc[i] = "Modal" + desc[i];
3395  }
3396 
3397  for (i = 6; i < 12; i++) {
3398  os.str(name);
3399  os.seekp(0, std::ios_base::end);
3400  os << ": " << sn_initial_dof[i/3] << xyz[i%3];
3401  desc[i] = os.str();
3402  }
3403 
3404  } else {
3405  if (i < 6 || i >= 12) {
3406  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3407  }
3408 
3409  os << ": " << sn_initial_dof[i/3] << xyz[i%3];
3410  desc[0] = os.str();
3411  }
3412 
3413 }
3414 
3415 std::ostream&
3416 ModalNode::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
3417 {
3418  if (bInitial == false) {
3419  integer iIndex = iGetFirstIndex();
3420 
3421  out
3422  << prefix << iIndex + 1 << "->" << iIndex + 3 << ": "
3423  "linear velocity definition [vx,vy,vz]" << std::endl
3424  << prefix << iIndex + 4 << "->" << iIndex + 6 << ": "
3425  "angular velocity definition [wx,wy,wz]" << std::endl;
3426  }
3427 
3428  StructNode::DescribeEq(out, prefix, bInitial);
3429 
3430  return out;
3431 }
3432 
3433 void
3434 ModalNode::DescribeEq(std::vector<std::string>& desc, bool bInitial, int i) const
3435 {
3436  if (bInitial || i == -1 || (i >= 6 && i < 12)) {
3437  int new_i = i;
3438  if (!bInitial && i != -1) {
3439  new_i = i - 6;
3440  }
3441  StructNode::DescribeEq(desc, bInitial, new_i);
3442 
3443  if (bInitial || (i >= 6 && i < 12)) {
3444  for (size_t ii = 0; ii < desc.size(); ii++) {
3445  desc[ii] = "Modal" + desc[ii];
3446  }
3447  return;
3448  }
3449  }
3450 
3451  if (i == -1) {
3452  desc.resize(12);
3453  for (int j = 0; j < 6; j++) {
3454  desc[6 + j] = "Modal" + desc[j];
3455  }
3456 
3457  } else {
3458  desc.resize(1);
3459  }
3460 
3461  std::ostringstream os;
3462  os << "ModalStructNode(" << GetLabel() << ")";
3463 
3464  const char **xeq = sn_modal_eq;
3465  if (bInitial) {
3466  xeq = sn_initial_eq;
3467  }
3468 
3469  if (i == -1) {
3470  std::string name = os.str();
3471 
3472  for (i = 0; i < 6; i++) {
3473  os.str(name);
3474  os.seekp(0, std::ios_base::end);
3475  os << ": " << xeq[i/3] << xyz[i%3];
3476  desc[i] = os.str();
3477  }
3478 
3479  } else {
3480  if (i < 0 || i >= 6) {
3481  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3482  }
3483 
3484  os << ": " << xeq[i/3] << xyz[i%3];
3485  desc[0] = os.str();
3486  }
3487 }
3488 
3489 /* Aggiorna dati in base alla soluzione */
3490 void
3491 ModalNode::Update(const VectorHandler& X, const VectorHandler& XP)
3492 {
3493  StructNode::Update(X, XP);
3494 
3495  integer iFirstIndex = iGetFirstIndex();
3496 
3497  /* aggiorno XPP e WP (servono solo a modal.cc) */
3498  XPPCurr = Vec3(XP, iFirstIndex + 7);
3499  WPCurr = Vec3(XP, iFirstIndex + 10);
3500 }
3501 
3502 void
3503 ModalNode::AfterConvergence(const VectorHandler& X,
3504  const VectorHandler& XP)
3505 {
3506  // override DynamicStructNode's function
3507  NO_OP;
3508 }
3509 
3510 /* ModalNode - end */
3511 
3512 
3513 /* DummyStructNode - begin */
3514 
3515 /* Costruttore definitivo */
3516 DummyStructNode::DummyStructNode(unsigned int uL,
3517  const DofOwner* pDO,
3518  const StructNode* pN,
3519  OrientationDescription ood,
3520  flag fOut)
3521 :
3522 StructDispNode(uL, pDO, ::Zero3, ::Zero3, 0, 0, 0., 0., ood, fOut),
3523 StructNode(uL, pDO, ::Zero3, ::Zero3x3, ::Zero3, ::Zero3, 0, 0, 0., 0., 0, ood, fOut),
3524 pNode(pN)
3525 {
3526  ASSERT(pNode != NULL);
3527 }
3528 
3529 
3530 /* Distruttore (per ora e' banale) */
3531 DummyStructNode::~DummyStructNode(void)
3532 {
3533  NO_OP;
3534 }
3535 
3536 
3537 /* Tipo di nodo strutturale */
3538 StructNode::Type
3539 DummyStructNode::GetStructNodeType(void) const
3540 {
3541  return StructNode::DUMMY;
3542 }
3543 
3544 StructDispNode::Type
3545 DummyStructNode::GetStructDispNodeType(void) const
3546 {
3547  return StructDispNode::UNKNOWN;
3548 }
3549 
3550 /* Restituisce il valore del dof iDof;
3551  * se differenziale, iOrder puo' essere = 1 per la derivata */
3552 const doublereal&
3553 DummyStructNode::dGetDofValue(int iDof, int iOrder) const
3554 {
3555  silent_cerr("DummyStructNode(" << GetLabel() << ") has no dofs"
3556  << std::endl);
3557  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3558 }
3559 
3560 
3561 /* Restituisce il valore del dof iDof al passo precedente;
3562  * se differenziale, iOrder puo' essere = 1 per la derivata */
3563 const doublereal&
3564 DummyStructNode::dGetDofValuePrev(int iDof, int iOrder) const
3565 {
3566  silent_cerr("DummyStructNode(" << GetLabel() << ") has no dofs"
3567  << std::endl);
3568  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3569 }
3570 
3571 
3572 /* Setta il valore del dof iDof a dValue;
3573  * se differenziale, iOrder puo' essere = 1 per la derivata */
3574 void
3575 DummyStructNode::SetDofValue(const doublereal& dValue,
3576  unsigned int iDof, unsigned int iOrder)
3577 {
3578  silent_cerr("DummyStructNode(" << GetLabel() << ") has no dofs"
3579  << std::endl);
3580  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3581 }
3582 
3583 
3584 /* Aggiorna dati durante l'iterazione fittizia iniziale */
3585 void
3586 DummyStructNode::DerivativesUpdate(const VectorHandler& X, const VectorHandler& XP)
3587 {
3588  /* posso farlo perche' in genere i dummy nodes si limitano
3589  * a copiare i valori di altri nodi, quindi non alterano
3590  * le variabili cinematiche */
3591  Update(X, XP);
3592 }
3593 
3594 
3595 /* Aggiorna dati in base alla soluzione durante l'assemblaggio iniziale */
3596 void
3597 DummyStructNode::InitialUpdate(const VectorHandler& /* X */ )
3598 {
3599  NO_OP;
3600 }
3601 
3602 
3603 /* Funzioni di inizializzazione, ereditate da DofOwnerOwner */
3604 void
3605 DummyStructNode::SetInitialValue(VectorHandler& /* X */ )
3606 {
3607  NO_OP;
3608 }
3609 
3610 
3611 void
3612 DummyStructNode::SetValue(DataManager *pDM,
3613  VectorHandler& X,
3614  VectorHandler& XP,
3615  SimulationEntity::Hints *ph)
3616 {
3617  Update(X, XP);
3618 }
3619 
3620 
3621 /* Elaborazione vettori e dati prima e dopo la predizione
3622  * per MultiStepIntegrator */
3623 void
3624 DummyStructNode::BeforePredict(VectorHandler& /* X */ ,
3625  VectorHandler& /* XP */ ,
3626  VectorHandler& /* XPrev */ ,
3627  VectorHandler& /* XPPrev */ ) const
3628 {
3629  NO_OP;
3630 }
3631 
3632 
3633 void
3634 DummyStructNode::AfterPredict(VectorHandler& X, VectorHandler& XP)
3635 {
3636  Update(X, XP);
3637 }
3638 
3639 bool
3640 DummyStructNode::ComputeAccelerations(bool b)
3641 {
3642  return const_cast<StructNode *>(pNode)->ComputeAccelerations(b);
3643 }
3644 
3645 /* DummyStructNode - end */
3646 
3647 
3648 /* OffsetDummyStructNode - begin */
3649 
3650 /* Costruttore definitivo */
3651 OffsetDummyStructNode::OffsetDummyStructNode(unsigned int uL,
3652  const DofOwner* pDO,
3653  const StructNode* pN,
3654  const Vec3& f,
3655  const Mat3x3& R,
3656  OrientationDescription ood,
3657  flag fOut)
3658 :
3659 StructDispNode(uL, pDO, ::Zero3, ::Zero3, 0, 0, 0., 0., ood, fOut),
3660 DummyStructNode(uL, pDO, pN, ood, fOut), f(f), R(R)
3661 {
3662  if (pNode->bOutputAccelerations()) {
3663  bOutputAccels = true;
3664  // const_cast<StructDispNode *>(this)->bOutputAccels = true;
3665 
3666  } else {
3667  bOutputAccels = false;
3668  // const_cast<StructDispNode *>(this)->bOutputAccels = false;
3669  }
3670 
3671  /* forzo la ricostruzione del nodo strutturale sottostante */
3672  Update_int();
3673 }
3674 
3675 
3676 /* Distruttore (per ora e' banale) */
3677 OffsetDummyStructNode::~OffsetDummyStructNode(void)
3678 {
3679  NO_OP;
3680 }
3681 
3682 
3683 /* update - interno */
3684 void
3685 OffsetDummyStructNode::Update_int(void)
3686 {
3687  Vec3 fCurr(pNode->GetRCurr()*f);
3688  XCurr = pNode->GetXCurr() + fCurr;
3689  RCurr = pNode->GetRCurr()*R;
3690  WCurr = pNode->GetWCurr();
3691  VCurr = pNode->GetVCurr() + WCurr.Cross(fCurr);
3692 
3693  if (bComputeAccelerations()) {
3694  WPCurr = pNode->GetWPCurr();
3695  XPPCurr = pNode->GetXPPCurr()
3696  + WCurr.Cross(WCurr.Cross(fCurr))
3697  + WPCurr.Cross(fCurr);
3698  }
3699 }
3700 
3701 
3702 /* Tipo di nodo dummy */
3703 DummyStructNode::Type
3704 OffsetDummyStructNode::GetDummyType(void) const
3705 {
3706  return DummyStructNode::OFFSET;
3707 }
3708 
3709 
3710 /* Aggiorna dati in base alla soluzione */
3711 void
3712 OffsetDummyStructNode::Update(const VectorHandler& /* X */ ,
3713  const VectorHandler& /* XP */ )
3714 {
3715  Update_int();
3716 }
3717 
3718 /* OffsetDummyStructNode - end */
3719 
3720 
3721 /* RelFrameDummyStructNode - begin */
3722 
3723 /* Costruttore definitivo */
3724 RelFrameDummyStructNode::RelFrameDummyStructNode(unsigned int uL,
3725  const DofOwner* pDO,
3726  const StructNode* pN,
3727  const StructNode* pNR,
3728  const Vec3& fh,
3729  const Mat3x3& Rh,
3730  OrientationDescription ood,
3731  flag fOut)
3732 :
3733 StructDispNode(uL, pDO, ::Zero3, ::Zero3, 0, 0, 0., 0., ood, fOut),
3734 DummyStructNode(uL, pDO, pN, ood, fOut),
3735 pNodeRef(pNR),
3736 RhT(Rh.Transpose()),
3737 fhT(RhT*fh)
3738 {
3739  ASSERT(pNodeRef != NULL);
3740 
3741  /*
3742  * Note: Rh is transposed from the beginning because it is
3743  * never used directly;
3744  * fh is premultiplied by Rh.Transpose() for the same reason
3745  *
3746  * Formulas:
3747  *
3748  * R = RhT * RrT * Rn
3749  * X = RhT * RrT * (Xn - Xr)
3750  * W = RhT * RrT * (Wn - Wr)
3751  * V = RhT * RrT * (Vn - Vr - Wr x (Xn - Xr))
3752  *
3753  * by defining
3754  *
3755  * Rn = Rr * Rh * R
3756  * Xn = Xr + Rr * (fh + Rh * X)
3757  *
3758  * and differentiating with respect to time
3759  */
3760 
3761  if (pNode->bOutputAccelerations()
3762  && pNodeRef->bOutputAccelerations())
3763  {
3764  bOutputAccels = true;
3765 
3766  } else {
3767  bOutputAccels = false;
3768  }
3769 
3770  /* forzo la ricostruzione del nodo strutturale sottostante */
3771  Update_int();
3772 }
3773 
3774 
3775 /* Distruttore (per ora e' banale) */
3776 RelFrameDummyStructNode::~RelFrameDummyStructNode(void)
3777 {
3778  NO_OP;
3779 }
3780 
3781 
3782 /* update - interno */
3783 void
3784 RelFrameDummyStructNode::Update_int(void)
3785 {
3786  Mat3x3 RT(RhT.MulMT(pNodeRef->GetRCurr()));
3787  Vec3 XRel(pNode->GetXCurr() - pNodeRef->GetXCurr());
3788 
3789  RCurr = RT*pNode->GetRCurr();
3790  XCurr = RT*XRel - fhT;
3791  WCurr = RT*(pNode->GetWCurr() - pNodeRef->GetWCurr());
3792 
3793  VCurr = RT*(pNode->GetVCurr()
3794  - pNodeRef->GetVCurr()
3795  - pNodeRef->GetWCurr().Cross(XRel));
3796 
3797  if (bComputeAccelerations()) {
3798  WPCurr = RT*(pNode->GetWPCurr() - pNodeRef->GetWPCurr()
3799  - pNodeRef->GetWCurr().Cross(pNode->GetWCurr()));
3800  XPPCurr = RT*(pNode->GetXPPCurr() - pNodeRef->GetXPPCurr()
3801  - pNodeRef->GetWPCurr().Cross(XRel)
3802  - (pNodeRef->GetWCurr()*2.).Cross(pNode->GetVCurr() - pNodeRef->GetVCurr())
3803  + pNodeRef->GetWCurr().Cross(pNodeRef->GetWCurr().Cross(XRel)));
3804  }
3805 }
3806 
3807 
3808 /* Tipo di nodo dummy */
3809 DummyStructNode::Type
3810 RelFrameDummyStructNode::GetDummyType(void) const
3811 {
3812  return DummyStructNode::RELATIVEFRAME;
3813 }
3814 
3815 
3816 /* Aggiorna dati in base alla soluzione */
3817 void
3818 RelFrameDummyStructNode::Update(const VectorHandler& /* X */ ,
3819  const VectorHandler& /* XP */ )
3820 {
3821  Update_int();
3822 }
3823 
3824 bool
3825 RelFrameDummyStructNode::ComputeAccelerations(bool b)
3826 {
3827  bool ok = true;
3828 
3829  if (!const_cast<StructNode *>(pNode)->ComputeAccelerations(b)) {
3830  ok = false;
3831  }
3832 
3833  if (!const_cast<StructNode *>(pNodeRef)->ComputeAccelerations(b)) {
3834  ok = false;
3835  }
3836 
3837  return ok;
3838 }
3839 
3840 /* RelFrameDummyStructNode - end */
3841 
3842 
3843 /* PivotRelFrameDummyStructNode - begin */
3844 
3845 /* Costruttore definitivo */
3846 PivotRelFrameDummyStructNode::PivotRelFrameDummyStructNode(unsigned int uL,
3847  const DofOwner* pDO,
3848  const StructNode* pN,
3849  const StructNode* pNR,
3850  const Vec3& fh,
3851  const Mat3x3& Rh,
3852  const StructNode* pNR2,
3853  const Vec3& fh2,
3854  const Mat3x3& Rh2,
3855  OrientationDescription ood,
3856  flag fOut)
3857 :
3858 StructDispNode(uL, pDO, ::Zero3, ::Zero3, 0, 0, 0., 0., ood, fOut),
3859 RelFrameDummyStructNode(uL, pDO, pN, pNR, fh, Rh, ood, fOut),
3860 pNodeRef2(pNR2), Rh2(Rh2), fh2(fh2)
3861 {
3862  ASSERT(pNodeRef2 != NULL);
3863 
3864  /*
3865  * Note: Rh is transposed from the beginning because it is
3866  * never used directly;
3867  * fh is premultiplied by Rh.Transpose() for the same reason
3868  *
3869  * Formulas:
3870  *
3871  * R = RhT * RrT * Rn
3872  * X = RhT * RrT * (Xn - Xr)
3873  * W = RhT * RrT * (Wn - Wr)
3874  * V = RhT * RrT * (Vn - Vr - Wr x (Xn - Xr))
3875  *
3876  * by defining
3877  *
3878  * Rn = Rr * Rh * R
3879  * Xn = Xr + Rr * (fh + Rh * X)
3880  *
3881  * and differentiating with respect to time
3882  */
3883 
3884  if (pNode->bOutputAccelerations()
3885  && pNodeRef->bOutputAccelerations()
3886  && pNodeRef2->bOutputAccelerations())
3887  {
3888  bOutputAccels = true;
3889 
3890  } else {
3891  bOutputAccels = false;
3892  }
3893 
3894  /* forzo la ricostruzione del nodo strutturale sottostante */
3895  Update_int();
3896 }
3897 
3898 
3899 /* Distruttore (per ora e' banale) */
3900 PivotRelFrameDummyStructNode::~PivotRelFrameDummyStructNode(void)
3901 {
3902  NO_OP;
3903 }
3904 
3905 
3906 /* update - interno */
3907 void
3908 PivotRelFrameDummyStructNode::Update_int(void)
3909 {
3910  RelFrameDummyStructNode::Update_int();
3911 
3912  Mat3x3 R2(pNodeRef2->GetRCurr()*Rh2);
3913 
3914  XCurr = pNodeRef2->GetRCurr()*(Rh2*XCurr + fh2);
3915 
3916  if (bComputeAccelerations()) {
3917  WPCurr = pNodeRef2->GetWPCurr()
3918  + pNodeRef2->GetWCurr().Cross(R2*WCurr)
3919  + R2*WPCurr;
3920  XPPCurr = pNodeRef2->GetXPPCurr()
3921  + pNodeRef2->GetWPCurr().Cross(XCurr)
3922  + pNodeRef2->GetWCurr().Cross(pNodeRef2->GetWCurr().Cross(XCurr))
3923  + (pNodeRef2->GetWCurr()*2.).Cross(R2*VCurr)
3924  + R2*XPPCurr;
3925  }
3926 
3927  WCurr = pNodeRef2->GetWCurr() + R2*WCurr;
3928  VCurr = pNodeRef2->GetVCurr()
3929  + pNodeRef2->GetWCurr().Cross(XCurr)
3930  + R2*VCurr;
3931  RCurr = R2*RCurr;
3932  XCurr += pNodeRef2->GetXCurr();
3933 }
3934 
3935 
3936 /* Tipo di nodo dummy */
3937 DummyStructNode::Type
3938 PivotRelFrameDummyStructNode::GetDummyType(void) const
3939 {
3940  return DummyStructNode::PIVOTRELATIVEFRAME;
3941 }
3942 
3943 
3944 /* Aggiorna dati in base alla soluzione */
3945 void
3946 PivotRelFrameDummyStructNode::Update(const VectorHandler& /* X */ ,
3947  const VectorHandler& /* XP */ )
3948 {
3949  Update_int();
3950 }
3951 
3952 bool
3953 PivotRelFrameDummyStructNode::ComputeAccelerations(bool b)
3954 {
3955  bool ok = true;
3956 
3957  if (!const_cast<StructNode *>(pNode)->ComputeAccelerations(b)) {
3958  ok = false;
3959  }
3960 
3961  if (!const_cast<StructNode *>(pNodeRef)->ComputeAccelerations(b)) {
3962  ok = false;
3963  }
3964 
3965  if (!const_cast<StructNode *>(pNodeRef2)->ComputeAccelerations(b)) {
3966  ok = false;
3967  }
3968 
3969  return ok;
3970 }
3971 
3972 /* RelFrameDummyStructNode - end */
3973 
3974 
3975 /* Legge un nodo strutturale */
3976 
3977 Node*
3978 ReadStructNode(DataManager* pDM,
3979  MBDynParser& HP,
3980  DofOwner* pDO,
3981  unsigned int uLabel)
3982 {
3983  DEBUGCOUT("Entering ReadStructNode(" << uLabel << ")" << std::endl);
3984 
3985  const char* sKeyWords[] = {
3986  "static" "displacement",
3987  "dynamic" "displacement",
3988  "static",
3989  "dynamic",
3990  "modal",
3991  "dummy",
3992  "offset",
3993  "relative" "frame",
3994  0
3995  };
3996 
3997  /* enum delle parole chiave */
3998  enum KeyWords {
3999  UNKNOWN = -1,
4000 
4001  STATIC_DISP = 0,
4002  DYNAMIC_DISP,
4003  STATIC,
4004  DYNAMIC,
4005  MODAL,
4006  DUMMY,
4007 
4008  OFFSET,
4009  RELATIVEFRAME,
4010 
4011  LASTKEYWORD
4012  };
4013 
4014  /* tabella delle parole chiave */
4015  KeyTable K(HP, sKeyWords);
4016 
4017  /* lettura dati specifici */
4018  KeyWords CurrType((KeyWords)HP.IsKeyWord());
4019 
4020  /*
4021  * explicit node type required; default is no longer "DYNAMIC"
4022  */
4023  if (CurrType == UNKNOWN) {
4024  silent_cerr("StructNode(" << uLabel << "): "
4025  "missing node type at line " << HP.GetLineData()
4026  << std::endl);
4027  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4028  }
4029 
4030 #ifdef DEBUG
4031  switch (CurrType) {
4032  case STATIC_DISP:
4033  std::cout << "Static structural displacement node" << std::endl;
4034  break;
4035  case DYNAMIC_DISP:
4036  std::cout << "Dynamic structural displacement node" << std::endl;
4037  break;
4038  case STATIC:
4039  std::cout << "Static structural node" << std::endl;
4040  break;
4041  case DYNAMIC:
4042  std::cout << "Dynamic structural node" << std::endl;
4043  break;
4044  case DUMMY:
4045  std::cout << "Dummy structural node" << std::endl;
4046  break;
4047  case MODAL:
4048  std::cout << "Modal node" << std::endl;
4049  break;
4050  default:
4051  std::cout << "Unknown structural node" << std::endl;
4052  break;
4053  }
4054 #endif /* DEBUG */
4055 
4056  StructDispNode* pNd = NULL;
4057  OrientationDescription od = UNKNOWN_ORIENTATION_DESCRIPTION;
4058  KeyWords DummyType = UNKNOWN;
4059  if (CurrType == DUMMY) {
4060  const StructNode* pNode = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
4061 
4062  DummyType = KeyWords(HP.GetWord());
4063  switch (DummyType) {
4064  case OFFSET: {
4065  ReferenceFrame RF(pNode);
4066  Vec3 f(HP.GetPosRel(RF));
4067  Mat3x3 R(HP.GetRotRel(RF));
4068 
4069  od = ReadOptionalOrientationDescription(pDM, HP);
4070 
4071  flag fOut = pDM->fReadOutput(HP, Node::STRUCTURAL);
4072  SAFENEWWITHCONSTRUCTOR(pNd,
4073  OffsetDummyStructNode,
4074  OffsetDummyStructNode(uLabel, pDO, pNode,
4075  f, R, od, fOut));
4076  } break;
4077 
4078  case RELATIVEFRAME: {
4079  const StructNode* pNodeRef = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
4080 
4081  ReferenceFrame RF(pNodeRef);
4082 
4083  Vec3 fh(Zero3);
4084  if (HP.IsKeyWord("position")) {
4085  fh = HP.GetPosRel(RF);
4086  }
4087 
4088  Mat3x3 Rh(Eye3);
4089  if (HP.IsKeyWord("orientation")) {
4090  Rh = HP.GetRotRel(RF);
4091 
4092  od = ReadOptionalOrientationDescription(pDM, HP);
4093  }
4094 
4095  const StructNode *pNodeRef2 = 0;
4096  Vec3 fh2(Zero3);
4097  Mat3x3 Rh2(Eye3);
4098  if (HP.IsKeyWord("pivot" "node")) {
4099  pNodeRef2 = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
4100 
4101  if (HP.IsKeyWord("position")) {
4102  fh2 = HP.GetPosRel(RF);
4103  }
4104 
4105  if (HP.IsKeyWord("orientation")) {
4106  Rh2 = HP.GetRotRel(RF);
4107  }
4108  }
4109 
4110  flag fOut = pDM->fReadOutput(HP, Node::STRUCTURAL);
4111 
4112  if (pNodeRef2) {
4113  SAFENEWWITHCONSTRUCTOR(pNd,
4114  PivotRelFrameDummyStructNode,
4115  PivotRelFrameDummyStructNode(uLabel, pDO,
4116  pNode, pNodeRef, fh, Rh,
4117  pNodeRef2, fh2, Rh2, od, fOut));
4118 
4119  } else {
4120  SAFENEWWITHCONSTRUCTOR(pNd,
4121  RelFrameDummyStructNode,
4122  RelFrameDummyStructNode(uLabel, pDO,
4123  pNode, pNodeRef, fh, Rh, od, fOut));
4124  }
4125  } break;
4126 
4127  default:
4128  silent_cerr("StructNode(" << uLabel << "): "
4129  "unknown dummy node type "
4130  "at line " << HP.GetLineData() << std::endl);
4131  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4132  }
4133 
4134  } else {
4135  bool bDisp(false);
4136 
4137  if (CurrType == STATIC_DISP || CurrType == DYNAMIC_DISP) {
4138  bDisp = true;
4139  }
4140 
4141  /* posizione (vettore di 3 elementi) */
4142  if (!HP.IsKeyWord("position")) {
4143  pedantic_cerr("StructNode(" << uLabel << "): "
4144  "missing keyword \"position\" at line "
4145  << HP.GetLineData() << std::endl);
4146  }
4147  Vec3 X0(HP.GetPosAbs(::AbsRefFrame));
4148  DEBUGCOUT("X0 =" << std::endl << X0 << std::endl);
4149 
4150  /* sistema di riferimento (trucco dei due vettori) */
4151  Mat3x3 R0;
4152  if (!bDisp) {
4153  if (!HP.IsKeyWord("orientation")) {
4154  pedantic_cerr("StructNode(" << uLabel << "): "
4155  "missing keyword \"orientation\" at line "
4156  << HP.GetLineData() << std::endl);
4157  }
4158  R0 = HP.GetRotAbs(::AbsRefFrame);
4159 
4160  od = ReadOptionalOrientationDescription(pDM, HP);
4161 
4162  DEBUGCOUT("R0 =" << std::endl << R0 << std::endl);
4163  }
4164 
4165  /* Velocita' iniziali (due vettori di 3 elementi, con la possibilita'
4166  * di usare "null" per porli uguali a zero) */
4167  if (!HP.IsKeyWord("velocity")) {
4168  pedantic_cerr("StructNode(" << uLabel << "): "
4169  "missing keyword \"velocity\" at line "
4170  << HP.GetLineData() << std::endl);
4171  }
4172  Vec3 XPrime0(HP.GetVelAbs(::AbsRefFrame, X0));
4173  DEBUGCOUT("Xprime0 =" << std::endl << XPrime0 << std::endl);
4174 
4175  Vec3 Omega0;
4176  if (!bDisp) {
4177  if (!HP.IsKeyWord("angular" "velocity")) {
4178  pedantic_cerr("StructNode(" << uLabel << "): "
4179  "missing keyword \"angular velocity\" at line "
4180  << HP.GetLineData() << std::endl);
4181  }
4182  Omega0 = HP.GetOmeAbs(::AbsRefFrame);
4183  DEBUGCOUT("Omega0 =" << std::endl << Omega0 << std::endl);
4184  }
4185 
4186  const StructNode *pRefNode = 0;
4187  if (HP.IsKeyWord("prediction" "node")) {
4188  switch (CurrType) {
4189  case STATIC_DISP:
4190  case DYNAMIC_DISP:
4191  case STATIC:
4192  case DYNAMIC:
4193  break;
4194 
4195  default:
4196  silent_cerr("StructNode(" << uLabel << "): "
4197  "prediction node allowed "
4198  "for static and dynamic nodes only, "
4199  "at line " << HP.GetLineData()
4200  << std::endl);
4201  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4202  }
4203  pRefNode = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
4204 
4205 #ifndef MBDYN_X_RELATIVE_PREDICTION
4206  silent_cerr("warning, relative prediction disabled; "
4207  "absolute prediction will be used" << std::endl);
4208 #endif /* ! MBDYN_X_RELATIVE_PREDICTION */
4209  }
4210 
4211  const RigidBodyKinematics *pRBK = pDM->pGetRBK();
4212 
4213  /* Rigidezza in assemblaggio diversa da quella di default
4214  * e flag di output */
4215  doublereal dPosStiff = pDM->dGetInitialPositionStiffness();
4216  doublereal dVelStiff = pDM->dGetInitialVelocityStiffness();
4217  bool bOmRot = pDM->bDoesOmegaRotate();
4218 
4219  if (HP.IsArg()) {
4220  if (HP.IsKeyWord("assembly")) {
4221  dPosStiff = HP.GetReal(dPosStiff);
4222  dVelStiff = HP.GetReal(dVelStiff);
4223 
4224  DEBUGCOUT("Initial position stiffness: " << dPosStiff << std::endl);
4225  DEBUGCOUT("Initial velocity stiffness: " << dVelStiff << std::endl);
4226 
4227  if (!bDisp) {
4228  bOmRot = HP.GetYesNoOrBool(bOmRot);
4229 
4230  DEBUGCOUT("Omega rotates? : " << (bOmRot ? "yes" : "no") << std::endl);
4231  }
4232  }
4233  }
4234 
4235  pDO->SetScale(pDM->dReadScale(HP, DofOwner::STRUCTURALNODE));
4236 
4237  flag fOut = pDM->fReadOutput(HP, Node::STRUCTURAL);
4238  switch (CurrType) {
4239  case DYNAMIC_DISP:
4240  case DYNAMIC:
4241  case MODAL:
4242  {
4243  bool bGotAccels(false);
4244  bool bAccels(pDM->bOutputAccelerations());
4245 
4246  bool bGotInertia(false);
4247  bool bInertia(false);
4248 
4249  while (HP.IsArg()) {
4250  if (HP.IsKeyWord("accelerations")) {
4251  if (bGotAccels) {
4252  silent_cerr("StructNode(" << uLabel << "): "
4253  "\"accelerations\" already set, "
4254  "repeated at line " << HP.GetLineData()
4255  << std::endl);
4256  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4257  }
4258  bGotAccels = true;
4259 
4260  if (HP.IsArg()) {
4261  if (HP.GetYesNoOrBool(false)) {
4262  bAccels = true;
4263 
4264  } else {
4265  bAccels = false;
4266  }
4267 
4268  } else {
4269  // deprecated
4270  silent_cout("StructNode(" << uLabel << "): "
4271  "warning, \"accelerations\" needs \"yes\" or \"no\" "
4272  "at line " << HP.GetLineData() << std::endl);
4273  bAccels = true;
4274  }
4275 
4276  } else if (HP.IsKeyWord("output" "inertia")) {
4277  if (bGotInertia) {
4278  silent_cerr("StructNode(" << uLabel << "): "
4279  "\"inertia\" already set, "
4280  "repeated at line " << HP.GetLineData()
4281  << std::endl);
4282  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4283  }
4284  bGotInertia = true;
4285 
4286  if (!HP.IsArg()) {
4287  silent_cerr("StructNode(" << uLabel << "): "
4288  "\"inertia\" needs \"yes\" or \"no\" "
4289  "at line " << HP.GetLineData() << std::endl);
4290  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4291  }
4292 
4293  if (HP.GetYesNoOrBool(true)) {
4294  bInertia = true;
4295 
4296  } else {
4297  bInertia = false;
4298  }
4299 
4300  } else {
4301  break;
4302  }
4303  }
4304 
4305  if (fOut) {
4306  // restore legacy behavior
4307  if (!bGotInertia) {
4308  bInertia = true;
4309  }
4310 
4311  if (bInertia) {
4312  fOut |= StructDispNode::OUTPUT_INERTIA;
4313  }
4314 
4315  if (bAccels) {
4316  fOut |= StructDispNode::OUTPUT_ACCELERATIONS;
4317  }
4318  }
4319  } break;
4320 
4321  default:
4322  break;
4323  }
4324 
4325  if (CurrType == DYNAMIC && pDM->bIsStaticModel()) {
4326  pedantic_cout("DynamicStructNode(" << uLabel << ") turned into static" << std::endl);
4327  CurrType = STATIC;
4328 
4329  } else if (CurrType == DYNAMIC_DISP && pDM->bIsStaticModel()) {
4330  pedantic_cout("DynamicStructDispNode(" << uLabel << ") turned into static" << std::endl);
4331  CurrType = STATIC_DISP;
4332  }
4333 
4334  /* Se non c'e' il punto e virgola finale */
4335  if (HP.IsArg()) {
4336  silent_cerr("ReadStructNode(" << uLabel << "): semicolon expected "
4337  "at line " << HP.GetLineData() << std::endl);
4338  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4339  }
4340 
4341  /* costruzione del nodo */
4342  switch (CurrType) {
4343  case STATIC_DISP:
4344  SAFENEWWITHCONSTRUCTOR(pNd, StaticStructDispNode,
4345  StaticStructDispNode(uLabel, pDO,
4346  X0,
4347  XPrime0,
4348  pRefNode,
4349  pRBK,
4350  dPosStiff, dVelStiff,
4351  od, fOut));
4352  break;
4353 
4354  case DYNAMIC_DISP:
4355  SAFENEWWITHCONSTRUCTOR(pNd, DynamicStructDispNode,
4356  DynamicStructDispNode(uLabel, pDO,
4357  X0,
4358  XPrime0,
4359  pRefNode,
4360  pRBK,
4361  dPosStiff, dVelStiff,
4362  od, fOut));
4363 
4364  /* Incrementa il numero di elementi automatici dei nodi dinamici */
4365  pDM->IncElemCount(Elem::AUTOMATICSTRUCTURAL);
4366  break;
4367 
4368  case STATIC:
4369  SAFENEWWITHCONSTRUCTOR(pNd, StaticStructNode,
4370  StaticStructNode(uLabel, pDO,
4371  X0, R0,
4372  XPrime0, Omega0,
4373  pRefNode,
4374  pRBK,
4375  dPosStiff, dVelStiff,
4376  bOmRot, od, fOut));
4377  break;
4378 
4379  case DYNAMIC:
4380  SAFENEWWITHCONSTRUCTOR(pNd, DynamicStructNode,
4381  DynamicStructNode(uLabel, pDO,
4382  X0, R0,
4383  XPrime0, Omega0,
4384  pRefNode,
4385  pRBK,
4386  dPosStiff, dVelStiff,
4387  bOmRot, od, fOut));
4388 
4389  /* Incrementa il numero di elementi automatici dei nodi dinamici */
4390  pDM->IncElemCount(Elem::AUTOMATICSTRUCTURAL);
4391  break;
4392 
4393  case MODAL:
4394  SAFENEWWITHCONSTRUCTOR(pNd, ModalNode,
4395  ModalNode(uLabel, pDO,
4396  X0, R0,
4397  XPrime0, Omega0,
4398  pRBK,
4399  dPosStiff, dVelStiff,
4400  bOmRot, od, fOut));
4401  break;
4402 
4403  default:
4404  ASSERT(false);
4405  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4406  }
4407  }
4408 
4409  std::ostream& out = pDM->GetLogFile();
4410 
4411  const char *description = "structural node: ";
4412 
4413  switch (CurrType){
4414  case DUMMY:
4415  switch (DummyType){
4416  case RELATIVEFRAME:
4417  description = "relative frame structural node: ";
4418  break;
4419 
4420  default:
4421  break;
4422  }
4423  break;
4424 
4425  default:
4426  break;
4427  }
4428 
4429  out << description << uLabel
4430  << " ", pNd->GetXCurr().Write(out, " ")
4431  << " ";
4432  const StructNode *pSN = dynamic_cast<const StructNode *>(pNd);
4433  if (pSN) {
4434  switch (od) {
4435  case EULER_123:
4436  out << "euler123 ",
4437  (MatR2EulerAngles123(pSN->GetRCurr())*dRaDegr).Write(out, " ");
4438  break;
4439 
4440  case EULER_313:
4441  out << "euler313 ",
4442  (MatR2EulerAngles313(pSN->GetRCurr())*dRaDegr).Write(out, " ");
4443  break;
4444 
4445  case EULER_321:
4446  out << "euler321 ",
4447  (MatR2EulerAngles321(pSN->GetRCurr())*dRaDegr).Write(out, " ");
4448  break;
4449 
4450  case ORIENTATION_VECTOR:
4451  out << "phi ",
4452  RotManip::VecRot(pSN->GetRCurr()).Write(out, " ");
4453  break;
4454 
4455  case ORIENTATION_MATRIX:
4456  out << "mat ",
4457  pSN->GetRCurr().Write(out, " ");
4458  break;
4459 
4460  default:
4461  /* impossible */
4462  break;
4463  }
4464 
4465  } else {
4466  switch (od) {
4467  case EULER_123:
4468  out << "euler123 ",
4469  ::Zero3.Write(out, " ");
4470  break;
4471 
4472  case EULER_313:
4473  out << "euler313 ",
4474  ::Zero3.Write(out, " ");
4475  break;
4476 
4477  case EULER_321:
4478  out << "euler321 ",
4479  ::Zero3.Write(out, " ");
4480  break;
4481 
4482  case ORIENTATION_VECTOR:
4483  out << "phi ",
4484  ::Zero3.Write(out, " ");
4485  break;
4486 
4487  case ORIENTATION_MATRIX:
4488  out << "mat ",
4489  ::Eye3.Write(out, " ");
4490  break;
4491 
4492  default:
4493  /* impossible */
4494  break;
4495  }
4496  }
4497 
4498  out << std::endl;
4499 
4500  ASSERT(pNd != NULL);
4501 
4502  return pNd;
4503 } /* End of ReadStructNode() */
4504 
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
sdn_initial_eq
static const char * sdn_initial_eq[]
Definition: strnode.cc:61
MBDynParser::RF
RFType RF
Definition: mbpar.h:166
CGR_Rot::MatGm1
const MatGm1_Manip MatGm1
Definition: matvec3.cc:647
StructNode::RPrev
Mat3x3 RPrev
Definition: strnode.h:728
StructNode::WPrev
Vec3 WPrev
Definition: strnode.h:745
DummyStructNode::DerivativesUpdate
virtual void DerivativesUpdate(const VectorHandler &X, const VectorHandler &XP)
Definition: strnode.cc:3586
StructNode::WPPrev
Vec3 WPPrev
Definition: strnode.h:750
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
grad::dGetValue
scalar_func_type dGetValue(scalar_func_type d)
Definition: gradient.h:2845
DynamicStructDispNode::iGetFirstRowIndex
virtual integer iGetFirstRowIndex(void) const
Definition: strnode.cc:1237
DynamicStructDispNode::GetStructDispNodeType
virtual StructDispNode::Type GetStructDispNodeType(void) const
Definition: strnode.cc:1101
StructDispNode::GetW
const Vec3 & GetW(void) const
Definition: strnode.cc:176
DynamicStructNode::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:3017
StaticStructDispNode::StaticStructDispNode
StaticStructDispNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Vec3 &V0, const StructNode *pRN, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, OrientationDescription od, flag fOut)
Definition: strnode.cc:1407
ModalNode::~ModalNode
virtual ~ModalNode(void)
Definition: strnode.cc:3326
DynamicStructDispNode::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:1171
Zero3
const Vec3 Zero3(0., 0., 0.)
Vec3::Cross
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
CGR_Rot::Param
const Param_Manip Param
Definition: matvec3.cc:644
DynamicStructNode::dGetDofValue
virtual const doublereal & dGetDofValue(int iDof, int iOrder=0) const
Definition: strnode.cc:3159
StructDispNode::dVelocityStiffness
doublereal dVelocityStiffness
Definition: strnode.h:115
Vec3::Write
std::ostream & Write(std::ostream &out, const char *sFill=" ") const
Definition: matvec3.cc:738
StructDispNode::GetStructDispNodeType
virtual StructDispNode::Type GetStructDispNodeType(void) const =0
flag
long int flag
Definition: mbdyn.h:43
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
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
DynamicStructDispNode::bComputeAccelerations
virtual bool bComputeAccelerations(void) const
Definition: strnode.h:588
DynamicStructDispNode::~DynamicStructDispNode
virtual ~DynamicStructDispNode(void)
Definition: strnode.cc:1095
StructNode::bOmegaRot
bool bOmegaRot
Definition: strnode.h:756
body.h
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
RelFrameDummyStructNode::bComputeAccelerations
virtual bool bComputeAccelerations(void) const
Definition: strnode.h:1855
Node
Definition: node.h:67
DataManager::dGetInitialVelocityStiffness
const doublereal & dGetInitialVelocityStiffness(void) const
Definition: dataman2.cc:117
AutomaticStructDispElem::GetGPCurr
virtual const Vec3 & GetGPCurr(void) const
Definition: autostr.h:57
StaticStructNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const
Definition: strnode.cc:3293
DynamicStructDispNode::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:1114
StructDispNode::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:200
OutputHandler::UseNetCDF
bool UseNetCDF(int out) const
Definition: output.cc:491
StructDispNode::GetWP
const Vec3 & GetWP(void) const
Definition: strnode.cc:188
StructNode::GetgCurr
virtual const Vec3 & GetgCurr(void) const
Definition: strnode.h:982
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
MBDynParser::GetRotAbs
Mat3x3 GetRotAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1857
DynamicStructDispNode::BeforePredict
virtual void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
Definition: strnode.cc:1288
StructDispNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: strnode.cc:145
StructNode::WPCurr
Vec3 WPCurr
Definition: strnode.h:749
StructDispNode::AfterPredict
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
Definition: strnode.cc:731
Vec3::Norm
doublereal Norm(void) const
Definition: matvec3.h:263
StructDispNode::InitialUpdate
virtual void InitialUpdate(const VectorHandler &X)
Definition: strnode.cc:594
sdn_initial_dof
static const char * sdn_initial_dof[]
Definition: strnode.cc:57
StructDispNode::XPrev
Vec3 XPrev
Definition: strnode.h:89
DynamicStructDispNode::SetDofValue
virtual void SetDofValue(const doublereal &dValue, unsigned int iDof, unsigned int iOrder=0)
Definition: strnode.cc:1361
StructDispNode::GetX
const Vec3 & GetX(void) const
Definition: strnode.cc:158
Write
std::ostream & Write(std::ostream &out, const FullMatrixHandler &m, const char *s, const char *s2)
Definition: fullmh.cc:376
DynamicStructDispNode::ComputeAccelerations
virtual bool ComputeAccelerations(bool b)
Definition: strnode.cc:1387
StructDispNode::SetInitialValue
virtual void SetInitialValue(VectorHandler &X)
Definition: strnode.cc:628
DynamicStructDispNode::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: strnode.cc:1278
DummyStructNode::~DummyStructNode
virtual ~DummyStructNode(void)
Definition: strnode.cc:3531
DynamicStructNode::GetGCurr
virtual const Vec3 & GetGCurr(void) const
Definition: strnode.cc:3103
PivotRelFrameDummyStructNode::~PivotRelFrameDummyStructNode
virtual ~PivotRelFrameDummyStructNode(void)
Definition: strnode.cc:3900
StructNode::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: strnode.cc:2332
DummyStructNode::SetInitialValue
virtual void SetInitialValue(VectorHandler &X)
Definition: strnode.cc:3605
MBDynParser::GetVelAbs
Vec3 GetVelAbs(const ReferenceFrame &rf, const Vec3 &x)
Definition: mbpar.cc:1482
OutputHandler::StrNodes
std::ostream & StrNodes(void) const
Definition: output.h:401
OffsetDummyStructNode::~OffsetDummyStructNode
virtual ~OffsetDummyStructNode(void)
Definition: strnode.cc:3677
Eye3
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
OrientationDescription
OrientationDescription
Definition: matvec3.h:1597
StructDispNode::~StructDispNode
virtual ~StructDispNode(void)
Definition: strnode.cc:138
DynamicStructNode::dGetDofValuePrev
virtual const doublereal & dGetDofValuePrev(int iDof, int iOrder=0) const
Definition: strnode.cc:3196
StaticStructNode::StaticStructNode
StaticStructNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Mat3x3 &R0, const Vec3 &V0, const Vec3 &W0, const StructNode *pRN, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, bool bOmRot, OrientationDescription ood, flag fOut)
Definition: strnode.cc:3261
DummyStructNode::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: strnode.cc:3612
DummyStructNode::dGetDofValue
virtual const doublereal & dGetDofValue(int iDof, int iOrder=0) const
Definition: strnode.cc:3553
grad::index_type
integer index_type
Definition: gradient.h:104
DynamicStructDispNode::dGetDofValuePrev
virtual const doublereal & dGetDofValuePrev(int iDof, int iOrder=0) const
Definition: strnode.cc:1336
StructDispNode::SetDofValue
virtual void SetDofValue(const doublereal &dValue, unsigned int iDof, unsigned int iOrder=0)
Definition: strnode.cc:422
PivotRelFrameDummyStructNode::ComputeAccelerations
virtual bool ComputeAccelerations(bool b)
Definition: strnode.cc:3953
DummyStructNode::AfterPredict
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
Definition: strnode.cc:3634
StructDispNode::iGetPrivDataIdx
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: strnode.cc:836
NO_OP
#define NO_OP
Definition: myassert.h:74
ReadOptionalOrientationDescription
OrientationDescription ReadOptionalOrientationDescription(DataManager *pDM, MBDynParser &HP)
Definition: dataman3.cc:2531
DummyStructNode::dGetDofValuePrev
virtual const doublereal & dGetDofValuePrev(int iDof, int iOrder=0) const
Definition: strnode.cc:3564
StructDispNode::GetXPP
const Vec3 & GetXPP(void) const
Definition: strnode.cc:182
StaticStructDispNode::GetStructDispNodeType
virtual StructDispNode::Type GetStructDispNodeType(void) const
Definition: strnode.cc:1429
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
ModalNode::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:3349
DynamicStructNode::GetWP
const Vec3 & GetWP(void) const
Definition: strnode.cc:2952
DynamicStructNode::DynamicStructNode
DynamicStructNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Mat3x3 &R0, const Vec3 &V0, const Vec3 &W0, const StructNode *pRN, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, bool bOmRot, OrientationDescription ood, flag fOut)
Definition: strnode.cc:2915
DataManager::bIsStaticModel
bool bIsStaticModel(void) const
Definition: dataman.h:482
StructNode::AfterPredict
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
Definition: strnode.cc:2492
StructNode::gPCurr
Vec3 gPCurr
Definition: strnode.h:735
ModalNode::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: strnode.cc:3503
DynamicStructNode::SetDofValue
virtual void SetDofValue(const doublereal &dValue, unsigned int iDof, unsigned int iOrder=0)
Definition: strnode.cc:3225
grad::fabs
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
AutomaticStructDispElem::AddInertia
virtual void AddInertia(const doublereal &dm)
Definition: autostr.cc:62
Mat3x3::GetVec
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
DynamicStructNode::~DynamicStructNode
virtual ~DynamicStructNode(void)
Definition: strnode.cc:2938
DummyStructNode::SetDofValue
virtual void SetDofValue(const doublereal &dValue, unsigned int iDof, unsigned int iOrder=0)
Definition: strnode.cc:3575
sdn_dof
static const char * sdn_dof[]
Definition: strnode.cc:49
StructDispNode::GetDofType
virtual DofOrder::Order GetDofType(unsigned int) const
Definition: strnode.cc:362
StructDispNode::XPPCurr
Vec3 XPPCurr
Definition: strnode.h:95
RotManip::VecRot
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
HighParser::GetYesNoOrBool
virtual bool GetYesNoOrBool(bool bDefval=false)
Definition: parser.cc:1038
StructNode::Output
virtual void Output(OutputHandler &OH) const
Definition: strnode.cc:1846
DynamicStructNode::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: strnode.cc:3129
StructDispNode::pRefRBK
const RigidBodyKinematics * pRefRBK
Definition: strnode.h:118
StaticStructNode::~StaticStructNode
virtual ~StaticStructNode(void)
Definition: strnode.cc:3285
Mat3x3::MulTV
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
DynamicStructNode::BeforePredict
virtual void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
Definition: strnode.cc:3143
func
void func(const T &u, const T &v, const T &w, doublereal e, T &f)
Definition: gradienttest.cc:309
StructNode::RCurr
Mat3x3 RCurr
Definition: strnode.h:730
AutomaticStructDispElem::ComputeAccelerations
virtual void ComputeAccelerations(Vec3 &XPP) const
Definition: autostr.cc:51
PivotRelFrameDummyStructNode::bComputeAccelerations
virtual bool bComputeAccelerations(void) const
Definition: strnode.h:1901
Zero3x3
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
DofOwnerOwner::pGetDofOwner
virtual const DofOwner * pGetDofOwner(void) const
Definition: dofown.h:113
ModalNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const
Definition: strnode.cc:3334
DataManager::dReadScale
doublereal dReadScale(MBDynParser &HP, enum DofOwner::Type t) const
Definition: dataman3.cc:1491
DynamicStructDispNode::GetBCurr
virtual const Vec3 & GetBCurr(void) const
Definition: strnode.cc:1253
DofOwner::SetScale
void SetScale(const doublereal &d)
Definition: dofown.cc:44
AutomaticStructDispElem::GetGCurr
virtual const Vec3 & GetGCurr(void) const
Definition: autostr.h:55
AbsRefFrame
const ReferenceFrame AbsRefFrame(0, Vec3(0., 0., 0), Mat3x3(1., 0., 0., 0., 1., 0., 0., 0., 1.), Vec3(0., 0., 0), Vec3(0., 0., 0), EULER_123)
DofOwner::dGetScale
doublereal dGetScale(void) const
Definition: dofown.cc:38
StructNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const =0
MatR2EulerAngles313
Vec3 MatR2EulerAngles313(const Mat3x3 &R)
Definition: matvec3.cc:927
StructDispNode::DerivativesUpdate
virtual void DerivativesUpdate(const VectorHandler &X, const VectorHandler &XP)
Definition: strnode.cc:582
AutomaticStructDispElem::GetBPCurr
virtual const Vec3 & GetBPCurr(void) const
Definition: autostr.h:56
DataManager::IncElemCount
void IncElemCount(Elem::Type type)
Definition: dataman2.cc:129
StructDispNode::pRefNode
const StructNode * pRefNode
Definition: strnode.h:98
sn_initial_eq
static const char * sn_initial_eq[]
Definition: strnode.cc:90
sdn_eq
static const char * sdn_eq[]
Definition: strnode.cc:53
PivotRelFrameDummyStructNode::pNodeRef2
const StructNode * pNodeRef2
Definition: strnode.h:1866
DummyStructNode::InitialUpdate
virtual void InitialUpdate(const VectorHandler &X)
Definition: strnode.cc:3597
MBDynParser::GetPosRel
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
StructNode::WCurr
Vec3 WCurr
Definition: strnode.h:747
StructNode::BeforePredict
virtual void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
Definition: strnode.cc:2389
StructDispNode::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: strnode.cc:639
MBDynParser::pDM
DataManager * pDM
Definition: mbpar.h:252
sn_modal_eq
static const char * sn_modal_eq[]
Definition: strnode.cc:78
StructDispNode::bOutputAccelerations
virtual bool bOutputAccelerations(void) const
Definition: strnode.h:439
false
Definition: mbdyn.h:76
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
StructNode::GetRPrev
virtual const Mat3x3 & GetRPrev(void) const
Definition: strnode.h:1000
PivotRelFrameDummyStructNode::PivotRelFrameDummyStructNode
PivotRelFrameDummyStructNode(unsigned int uL, const DofOwner *pDO, const StructNode *pNode, const StructNode *pNodeRef, const Vec3 &fh, const Mat3x3 &Rh, const StructNode *pNodeRef2, const Vec3 &fh2, const Mat3x3 &Rh2, OrientationDescription ood, flag fOut)
Definition: strnode.cc:3846
StructDispNode::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: strnode.cc:341
sn_dof
static const char * sn_dof[]
Definition: strnode.cc:66
StructDispNode::od
OrientationDescription od
Definition: strnode.h:111
RigidBodyKinematics::Update
virtual void Update(void)
Definition: rbk.cc:45
CGR_Rot::MatG
const MatG_Manip MatG
Definition: matvec3.cc:646
MBDynParser::GetOmeAbs
Vec3 GetOmeAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1551
OutputHandler::GetCurrentStep
long GetCurrentStep(void) const
Definition: output.h:116
StructDispNode::BeforePredict
virtual void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
Definition: strnode.cc:680
grad::FunctionCall
FunctionCall
Definition: gradient.h:1018
DummyStructNode::DummyStructNode
DummyStructNode(unsigned int uL, const DofOwner *pDO, const StructNode *pNode, OrientationDescription ood, flag fOut)
Definition: strnode.cc:3516
AutomaticStructDispElem::GetBCurr
virtual const Vec3 & GetBCurr(void) const
Definition: autostr.h:54
DynamicStructDispNode::pAutoStr
AutomaticStructDispElem * pAutoStr
Definition: strnode.h:484
grad::Gradient::SetValuePreserve
void SetValuePreserve(scalar_func_type dVal)
Definition: gradient.h:2511
DataManager::bDoesOmegaRotate
bool bDoesOmegaRotate(void) const
Definition: dataman2.cc:123
StructNode::RRef
Mat3x3 RRef
Definition: strnode.h:729
Node::Type
Type
Definition: node.h:71
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
MBDynParser::GetPosAbs
Vec3 GetPosAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1401
MatR2EulerAngles123
Vec3 MatR2EulerAngles123(const Mat3x3 &R)
Definition: matvec3.cc:893
StructNode::GetR
const Mat3x3 & GetR(void) const
Definition: strnode.cc:1481
true
Definition: mbdyn.h:77
RelFrameDummyStructNode::~RelFrameDummyStructNode
virtual ~RelFrameDummyStructNode(void)
Definition: strnode.cc:3776
DummyStructNode::pNode
const StructNode * pNode
Definition: strnode.h:1660
RelFrameDummyStructNode::pNodeRef
const StructNode * pNodeRef
Definition: strnode.h:1823
StructDispNode::GetV
const Vec3 & GetV(void) const
Definition: strnode.cc:170
DummyStructNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const
Definition: strnode.cc:3539
xyz
static const char xyz[]
Definition: strnode.cc:47
StructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const =0
StructNode::GetWCurr
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
StructDispNode::XPPPrev
Vec3 XPPPrev
Definition: strnode.h:96
dRaDegr
const doublereal dRaDegr
Definition: matvec3.cc:884
StructNode::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:1562
StructNode::DerivativesUpdate
virtual void DerivativesUpdate(const VectorHandler &X, const VectorHandler &XP)
Definition: strnode.cc:2240
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
StructDispNode::bOutputAccels
bool bOutputAccels
Definition: strnode.h:122
StructNode::GetWPCurr
virtual const Vec3 & GetWPCurr(void) const
Definition: strnode.h:1042
OutputHandler::IsOpen
bool IsOpen(int out) const
Definition: output.cc:395
StructDispNode::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP, const VectorHandler &XPP)
Definition: strnode.cc:789
CGR_Rot::MatR
const MatR_Manip MatR
Definition: matvec3.cc:645
DummyStructNode::GetStructDispNodeType
virtual StructDispNode::Type GetStructDispNodeType(void) const
Definition: strnode.cc:3545
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
OffsetDummyStructNode::OffsetDummyStructNode
OffsetDummyStructNode(unsigned int uL, const DofOwner *pDO, const StructNode *pNode, const Vec3 &f, const Mat3x3 &R, OrientationDescription ood, flag fOut)
Definition: strnode.cc:3651
StructDispNode::pGetRBK
const RigidBodyKinematics * pGetRBK(void) const
Definition: strnode.cc:152
KeyWords
KeyWords
Definition: dataman4.cc:94
Mat3x3::MulTM
Mat3x3 MulTM(const Mat3x3 &m) const
Definition: matvec3.cc:500
StructDispNode::dGetPrivData
virtual doublereal dGetPrivData(unsigned int i) const
Definition: strnode.cc:980
sn_initial_dof
static const char * sn_initial_dof[]
Definition: strnode.cc:84
StructNode::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP, const VectorHandler &XPP)
Definition: strnode.cc:2592
StructDispNode::dPositionStiffness
doublereal dPositionStiffness
Definition: strnode.h:114
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
StructNode::iGetPrivDataIdx
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: strnode.cc:2662
grad::Gradient::DerivativeResizeReset
void DerivativeResizeReset(LocalDofMap *pMap, index_type iStartGlobal, index_type iEndGlobal, MapVectorBase::GlobalScope s, scalar_deriv_type dVal)
Definition: gradient.h:2538
StructNode::gPRef
Vec3 gPRef
Definition: strnode.h:734
StructNode::InitialUpdate
virtual void InitialUpdate(const VectorHandler &X)
Definition: strnode.cc:2258
StructNode::WRef
Vec3 WRef
Definition: strnode.h:746
StructNode::~StructNode
virtual ~StructNode(void)
Definition: strnode.cc:1475
RelFrameDummyStructNode::GetDummyType
virtual DummyStructNode::Type GetDummyType(void) const
Definition: strnode.cc:3810
StructDispNode::ComputeAccelerations
virtual bool ComputeAccelerations(bool b)
Definition: strnode.cc:194
StructNode::OutputPrepare
virtual void OutputPrepare(OutputHandler &OH)
Definition: strnode.cc:1779
DynamicStructDispNode::GetBPCurr
virtual const Vec3 & GetBPCurr(void) const
Definition: strnode.cc:1259
StructNode::gRef
Vec3 gRef
Definition: strnode.h:732
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
SimulationEntity::iGetNumDof
virtual unsigned int iGetNumDof(void) const =0
DynamicStructDispNode::GetXPP
const Vec3 & GetXPP(void) const
Definition: strnode.cc:1108
DummyStructNode::ComputeAccelerations
virtual bool ComputeAccelerations(bool b)
Definition: strnode.cc:3640
DynamicStructNode::GetGPCurr
virtual const Vec3 & GetGPCurr(void) const
Definition: strnode.cc:3109
dS
const doublereal dS
Definition: beamslider.cc:71
Mat3x3::pGetMat
const doublereal * pGetMat(void) const
Definition: matvec3.h:743
VectorHandler::Put
virtual void Put(integer iRow, const Vec3 &v)
Definition: vh.cc:93
StructDispNode::dGetDofValue
virtual const doublereal & dGetDofValue(int iDof, int iOrder=0) const
Definition: strnode.cc:372
DataManager::pGetRBK
const RigidBodyKinematics * pGetRBK(void) const
Definition: dataman.h:485
RelFrameDummyStructNode::ComputeAccelerations
virtual bool ComputeAccelerations(bool b)
Definition: strnode.cc:3825
sn_eq
static const char * sn_eq[]
Definition: strnode.cc:72
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
StructNode::SetDofValue
virtual void SetDofValue(const doublereal &dValue, unsigned int iDof, unsigned int iOrder=0)
Definition: strnode.cc:1746
StructDispNode::VCurr
Vec3 VCurr
Definition: strnode.h:93
StructNode::dGetPrivData
virtual doublereal dGetPrivData(unsigned int i) const
Definition: strnode.cc:2806
StaticStructDispNode::~StaticStructDispNode
virtual ~StaticStructDispNode(void)
Definition: strnode.cc:1423
StructNode::dGetDofValuePrev
virtual const doublereal & dGetDofValuePrev(int iDof, int iOrder=0) const
Definition: strnode.cc:1713
StructDispNode::OutputPrepare
virtual void OutputPrepare(OutputHandler &OH)
Definition: strnode.cc:445
StructNode::iGetNumPrivData
virtual unsigned int iGetNumPrivData(void) const
Definition: strnode.cc:2630
DataManager::dGetInitialPositionStiffness
const doublereal & dGetInitialPositionStiffness(void) const
Definition: dataman2.cc:111
HighParser::GetWord
virtual int GetWord(void)
Definition: parser.cc:1083
ReadStructNode
Node * ReadStructNode(DataManager *pDM, MBDynParser &HP, DofOwner *pDO, unsigned int uLabel)
Definition: strnode.cc:3978
Vec3::pGetVec
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
DataManager::GetLogFile
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
ModalNode::iGetFirstRowIndex
virtual integer iGetFirstRowIndex(void) const
Definition: strnode.cc:3343
DummyStructNode::BeforePredict
virtual void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
Definition: strnode.cc:3624
MatR2EulerAngles321
Vec3 MatR2EulerAngles321(const Mat3x3 &R)
Definition: matvec3.cc:941
OffsetDummyStructNode::Update_int
void Update_int(void)
Definition: strnode.cc:3685
StructNode::GetgPCurr
virtual const Vec3 & GetgPCurr(void) const
Definition: strnode.h:994
DynamicStructDispNode::DynamicStructDispNode
DynamicStructDispNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Vec3 &V0, const StructNode *pRN, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, OrientationDescription od, flag fOut)
Definition: strnode.cc:1076
StructDispNode::XCurr
Vec3 XCurr
Definition: strnode.h:90
StructDispNode::dGetDofValuePrev
virtual const doublereal & dGetDofValuePrev(int iDof, int iOrder=0) const
Definition: strnode.cc:397
StructDispNode::GetVCurr
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
DynamicStructNode::iGetFirstRowIndex
virtual integer iGetFirstRowIndex(void) const
Definition: strnode.cc:3085
ModalNode::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:3416
StructNode::SetInitialValue
virtual void SetInitialValue(VectorHandler &X)
Definition: strnode.cc:2319
ModalNode::ModalNode
ModalNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Mat3x3 &R0, const Vec3 &V0, const Vec3 &W0, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, bool bOmRot, OrientationDescription ood, flag fOut)
Definition: strnode.cc:3303
StructNode::StructNode
StructNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Mat3x3 &R0, const Vec3 &V0, const Vec3 &W0, const StructNode *pRN, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, bool bOmRot, OrientationDescription ood, flag fOut)
Definition: strnode.cc:1439
OffsetDummyStructNode::GetDummyType
virtual DummyStructNode::Type GetDummyType(void) const
Definition: strnode.cc:3704
Mat3x3::MulMT
Mat3x3 MulMT(const Mat3x3 &m) const
Definition: matvec3.cc:444
StructDispNode::GetXPPCurr
virtual const Vec3 & GetXPPCurr(void) const
Definition: strnode.h:334
StructDispNode::StructDispNode
StructDispNode(unsigned int uL, const DofOwner *pDO, const Vec3 &X0, const Vec3 &V0, const StructNode *pRN, const RigidBodyKinematics *pRBK, doublereal dPosStiff, doublereal dVelStiff, OrientationDescription od, flag fOut)
Definition: strnode.cc:102
DynamicStructDispNode::dGetDofValue
virtual const doublereal & dGetDofValue(int iDof, int iOrder=0) const
Definition: strnode.cc:1303
StructNode::dGetDofValue
virtual const doublereal & dGetDofValue(int iDof, int iOrder=0) const
Definition: strnode.cc:1680
ToBeOutput::SetOutputFlag
virtual void SetOutputFlag(flag f=flag(1))
Definition: output.cc:896
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
StructNode::GetW
const Vec3 & GetW(void) const
Definition: strnode.cc:1487
StructDispNode::VPrev
Vec3 VPrev
Definition: strnode.h:92
DynamicStructNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const
Definition: strnode.h:1373
StructNode::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: strnode.cc:1652
doublereal
double doublereal
Definition: colamd.c:52
StructNode::GetWP
const Vec3 & GetWP(void) const
Definition: strnode.cc:1493
StructNode::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:1499
integer
long int integer
Definition: colamd.c:51
DynamicStructDispNode::AddInertia
virtual void AddInertia(const doublereal &dm) const
Definition: strnode.cc:1244
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
Mat3x3::Write
std::ostream & Write(std::ostream &out, const char *sFill=" ", const char *sFill2=NULL) const
Definition: matvec3.cc:722
MatR2EulerParams
void MatR2EulerParams(const Mat3x3 &R, doublereal &e0, Vec3 &e)
Definition: matvec3.cc:954
DynamicStructNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const
Definition: strnode.cc:2946
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
RelFrameDummyStructNode::RhT
const Mat3x3 RhT
Definition: strnode.h:1824
GRADIENT_ASSERT
#define GRADIENT_ASSERT(expr)
Definition: gradient.h:74
DynamicStructDispNode::SetOutputFlag
virtual void SetOutputFlag(flag f=flag(1))
Definition: strnode.cc:1394
StructDispNode::iGetNumPrivData
virtual unsigned int iGetNumPrivData(void) const
Definition: strnode.cc:804
DataManager::ReadNode
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
RelFrameDummyStructNode::RelFrameDummyStructNode
RelFrameDummyStructNode(unsigned int uL, const DofOwner *pDO, const StructNode *pNode, const StructNode *pNodeRef, const Vec3 &fh, const Mat3x3 &Rh, OrientationDescription ood, flag fOut)
Definition: strnode.cc:3724
StructDispNode::GetR
const Mat3x3 & GetR(void) const
Definition: strnode.cc:164
PivotRelFrameDummyStructNode::GetDummyType
virtual DummyStructNode::Type GetDummyType(void) const
Definition: strnode.cc:3938
StructDispNode::bComputeAccelerations
virtual bool bComputeAccelerations(void) const
Definition: strnode.h:433
R
Mat3x3 R
Definition: mbdyn_util_oct.cc:74
StructDispNode::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:259
DynamicStructNode::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:2958
OutputHandler::UseText
bool UseText(int out) const
Definition: output.cc:446
DynamicStructNode::AddInertia
virtual void AddInertia(const doublereal &dm, const Vec3 &dS, const Mat3x3 &dJ) const
Definition: strnode.cc:3092
DataManager::bOutputAccelerations
bool bOutputAccelerations(void) const
Definition: dataman.cc:878
StructNode::gCurr
Vec3 gCurr
Definition: strnode.h:733
DynamicStructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const
Definition: strnode.h:602
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056