MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
beam.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/beam.cc,v 1.110 2017/05/12 17:29:26 morandini Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /*
33  * Trave a volumi finiti, con predisposizione per forze di inerzia consistenti
34  * e legame cositutivo piezoelettico
35  */
36 
37 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
38 
39 #include <limits>
40 #include <cfloat>
41 #include <limits>
42 
43 #include "dataman.h"
44 #include "constltp.h"
45 #include "shapefnc.h"
46 #include "beam.h"
47 #include "pzbeam.h"
48 #include "Rot.hh"
49 
50 /*
51  * Nota: non e' ancora stato implementato il contributo
52  * della ViscoElasticBeam all'assemblaggio iniziale
53  */
54 
55 /*
56  * Nota: la parte viscoelastica va rivista in accordo con la piu'
57  * recente formulazione delle derivate delle deformazioni nel sistema
58  * materiale
59  */
60 
61 /* Beam - begin */
62 
63 /* Costruttore normale */
64 Beam::Beam(unsigned int uL,
65  const StructNode* pN1,
66  const StructNode* pN2,
67  const StructNode* pN3,
68  const Vec3& F1,
69  const Vec3& F2,
70  const Vec3& F3,
71  const Mat3x3& R1,
72  const Mat3x3& R2,
73  const Mat3x3& R3,
74  const Mat3x3& r_I,
75  const Mat3x3& rII,
76  const ConstitutiveLaw6D* pD_I,
77  const ConstitutiveLaw6D* pDII,
78  OrientationDescription ood,
79  flag fOut)
80 : Elem(uL, fOut),
81 ElemGravityOwner(uL, fOut),
82 InitialAssemblyElem(uL, fOut),
83 od(ood),
84 f(),
85 RNode(),
86 bConsistentInertia(false),
87 dMass_I(0.),
88 S0_I(Zero3),
89 J0_I(Zero3x3),
90 dMassII(0.),
91 S0II(Zero3),
92 J0II(Zero3x3),
93 bFirstRes(false)
94 {
95  ASSERT(pN1 != NULL);
96  ASSERT(pN1->GetNodeType() == Node::STRUCTURAL);
97  ASSERT(pN2 != NULL);
98  ASSERT(pN2->GetNodeType() == Node::STRUCTURAL);
99  ASSERT(pN3 != NULL);
100  ASSERT(pN3->GetNodeType() == Node::STRUCTURAL);
101 
102  pNode[NODE1] = pN1;
103  pNode[NODE2] = pN2;
104  pNode[NODE3] = pN3;
105 
106  const_cast<Vec3&>(f[NODE1]) = F1;
107  const_cast<Vec3&>(f[NODE2]) = F2;
108  const_cast<Vec3&>(f[NODE3]) = F3;
109  const_cast<Mat3x3&>(RNode[NODE1]) = R1;
110  const_cast<Mat3x3&>(RNode[NODE2]) = R2;
111  const_cast<Mat3x3&>(RNode[NODE3]) = R3;
112  RRef[S_I] = R[S_I] = (Mat3x3&)r_I;
113  RRef[SII] = R[SII] = (Mat3x3&)rII;
114 
115  pD[S_I] = 0;
116  SAFENEWWITHCONSTRUCTOR(pD[S_I],
117  ConstitutiveLaw6DOwner,
118  ConstitutiveLaw6DOwner(pD_I));
119  pD[SII] = 0;
120  SAFENEWWITHCONSTRUCTOR(pD[SII],
121  ConstitutiveLaw6DOwner,
122  ConstitutiveLaw6DOwner(pDII));
123 
124  Init();
125 }
126 
127 
128 /* Costruttore per la trave con forze d'inerzia consistenti */
129 Beam::Beam(unsigned int uL,
130  const StructNode* pN1,
131  const StructNode* pN2,
132  const StructNode* pN3,
133  const Vec3& F1,
134  const Vec3& F2,
135  const Vec3& F3,
136  const Mat3x3& R1,
137  const Mat3x3& R2,
138  const Mat3x3& R3,
139  const Mat3x3& r_I,
140  const Mat3x3& rII,
141  const ConstitutiveLaw6D* pD_I,
142  const ConstitutiveLaw6D* pDII,
143  doublereal dM_I,
144  const Vec3& s0_I,
145  const Mat3x3& j0_I,
146  doublereal dMII,
147  const Vec3& s0II,
148  const Mat3x3& j0II,
149  OrientationDescription ood,
150  flag fOut)
151 : Elem(uL, fOut),
152 ElemGravityOwner(uL, fOut),
153 InitialAssemblyElem(uL, fOut),
154 od(ood),
155 f(),
156 RNode(),
157 bConsistentInertia(true),
158 dMass_I(dM_I),
159 S0_I(s0_I),
160 J0_I(j0_I),
161 dMassII(dMII),
162 S0II(s0II),
163 J0II(j0II),
164 bFirstRes(false)
165 {
166  pNode[NODE1] = pN1;
167  pNode[NODE2] = pN2;
168  pNode[NODE3] = pN3;
169  const_cast<Vec3&>(f[NODE1]) = F1;
170  const_cast<Vec3&>(f[NODE2]) = F2;
171  const_cast<Vec3&>(f[NODE3]) = F3;
172  const_cast<Mat3x3&>(RNode[NODE1]) = R1;
173  const_cast<Mat3x3&>(RNode[NODE2]) = R2;
174  const_cast<Mat3x3&>(RNode[NODE3]) = R3;
175  RPrev[S_I] = RRef[S_I] = R[S_I] = const_cast<Mat3x3&>(r_I);
176  RPrev[SII] = RRef[SII] = R[SII] = const_cast<Mat3x3&>(rII);
177 
178  pD[S_I] = 0;
179  SAFENEWWITHCONSTRUCTOR(pD[S_I],
180  ConstitutiveLaw6DOwner,
181  ConstitutiveLaw6DOwner(pD_I));
182  pD[SII] = 0;
183  SAFENEWWITHCONSTRUCTOR(pD[SII],
184  ConstitutiveLaw6DOwner,
185  ConstitutiveLaw6DOwner(pDII));
186 
187  Init();
188 }
189 
190 void
191 Beam::Init(void)
192 {
193  for (unsigned i = 0; i < NUMSEZ; i++) {
194 #ifdef USE_NETCDF
195  Var_X[i] = 0;
196  Var_Phi[i] = 0;
197  Var_F[i] = 0;
198  Var_M[i] = 0;
199  Var_Nu[i] = 0;
200  Var_K[i] = 0;
201  Var_NuP[i] = 0;
202  Var_KP[i] = 0;
203 #endif /* USE_NETCDF */
204 
205  Omega[i] = Zero3;
206  Az[i] = Zero6;
207  AzRef[i] = Zero6;
208  AzLoc[i] = Zero6;
209  DefLoc[i] = Zero6;
210  DefLocRef[i] = Zero6;
211  DefLocPrev[i] = Zero6;
212  p[i] = Zero3;
213  g[i] = Zero3;
214  L0[i] = Zero3;
215  L[i] = Zero3;
216  LRef[i] = Zero3;
217  }
218 
219  DsDxi();
220 
221  Vec3 xTmp[NUMNODES];
222 
223  for (unsigned int i = 0; i < NUMNODES; i++) {
224  xTmp[i] = pNode[i]->GetXCurr()+pNode[i]->GetRCurr()*f[i];
225  }
226 
227  /* Aggiorna le grandezze della trave nei punti di valutazione */
228  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
229  p[iSez] = InterpState(xTmp[NODE1],
230  xTmp[NODE2],
231  xTmp[NODE3],
232  Beam::Section(iSez));
233  }
234 }
235 
236 Beam::~Beam(void)
237 {
238  ASSERT(pD[S_I] != 0);
239  if (pD[S_I] != 0) {
240  SAFEDELETE(pD[S_I]);
241  }
242 
243  ASSERT(pD[SII] != 0);
244  if (pD[SII] != 0) {
245  SAFEDELETE(pD[SII]);
246  }
247 }
248 
249 static const unsigned int iNumPrivData =
250  +3 // 0 ( 1 -> 3) - "e" strain
251  +3 // 3 ( 4 -> 6) - "k" curvature
252  +3 // 6 ( 7 -> 9) - "F" force
253  +3 // 9 (10 -> 12) - "M" moment
254  +3 // 12 (13 -> 15) - "X" position
255  +3 // 15 (16 -> 18) - "Phi" orientation vector
256  +3 // 18 (19 -> 21) - "Omega" angular velocity
257  +3 // 21 (22 -> 24) - "eP" strain rate
258  +3 // 24 (25 -> 27) - "kP" curvature rate
259  ;
260 
261 /* Accesso ai dati privati */
262 unsigned int
263 Beam::iGetNumPrivData(void) const
264 {
265  return 2*iNumPrivData;
266 }
267 
268 unsigned int
269 Beam::iGetPrivDataIdx_int(const char *s, ConstLawType::Type type)
270 {
271  ASSERT(s != NULL);
272 
273  unsigned int idx = 0;
274 
275  switch (s[0]) {
276  case 'k':
277  idx += 3;
278  case 'e':
279  if (s[1] == 'P') {
280  if (type != ConstLawType::VISCOUS) {
281  return 0;
282  }
283  s++;
284  idx += 21;
285  break;
286  }
287  break;
288 
289  case 'F':
290  idx += 6;
291  break;
292 
293  case 'M':
294  idx += 9;
295  break;
296 
297  case 'X':
298  idx += 12;
299  break;
300 
301  case 'P':
302  if (strncasecmp(s, "Phi", STRLENOF("Phi")) != 0) {
303  return 0;
304  }
305  s += STRLENOF("Phi") - 1;
306  idx += 15;
307  break;
308 
309  case 'O':
310  if (strncasecmp(s, "Omega", STRLENOF("Omega")) != 0) {
311  return 0;
312  }
313  s += STRLENOF("Omega") - 1;
314  idx += 18;
315  break;
316 
317  default:
318  return 0;
319  }
320 
321  switch (s[1]) {
322  case 'x':
323  idx += 1;
324  break;
325 
326  case 'y':
327  idx += 2;
328  break;
329 
330  case 'z':
331  idx += 3;
332  break;
333 
334  default:
335  return 0;
336  }
337 
338  if (s[2] != '\0') {
339  return 0;
340  }
341 
342  return idx;
343 }
344 
345 unsigned int
346 Beam::iGetPrivDataIdx(const char *s) const
347 {
348  ASSERT(s != NULL);
349 
350  unsigned int p_idx = 0;
351 
352  if (strncmp(s, "pI", STRLENOF("pI")) != 0) {
353  return 0;
354  }
355  s += STRLENOF("pI");
356 
357  if (s[0] == 'I') {
358  p_idx += iNumPrivData;
359  s++;
360  }
361 
362  if (s[0] != '.') {
363  return 0;
364  }
365  s++;
366 
367  ConstLawType::Type type = ConstLawType::ELASTIC;
368  if (dynamic_cast<const ViscoElasticBeam *>(this)) {
369  type = ConstLawType::VISCOUS;
370  }
371 
372  unsigned int idx = iGetPrivDataIdx_int(s, type);
373  if (idx != 0) {
374  idx += p_idx;
375  }
376 
377  return idx;
378 }
379 
380 doublereal
381 Beam::dGetPrivData(unsigned int i) const
382 {
383  ASSERT(i > 0 && i <= iGetNumPrivData());
384 
385  int sez = (i - 1)/iNumPrivData;
386  int idx = (i - 1)%iNumPrivData + 1;
387 
388  switch (idx) {
389  case 1:
390  case 2:
391  case 3:
392  case 4:
393  case 5:
394  case 6:
395 
396  return DefLoc[sez].dGet(idx);
397 
398  case 7:
399  case 8:
400  case 9:
401  case 10:
402  case 11:
403  case 12:
404  return AzLoc[sez].dGet(idx - 6);
405 
406  case 13:
407  case 14:
408  case 15:
409  return p[sez].dGet(idx - 12);
410 
411  case 16:
412  case 17:
413  case 18:
414  return RotManip::VecRot(R[sez]).dGet(idx - 15);
415 
416  case 19:
417  case 20:
418  case 21:
419  return Omega[sez].dGet(idx - 18);
420 
421  default:
422  silent_cerr("Beam(" << GetLabel() << "): "
423  "illegal private data " << i << std::endl);
424  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
425  }
426 }
427 
428 Vec3
429 Beam::InterpState(const Vec3& v1,
430  const Vec3& v2,
431  const Vec3& v3,
432  enum Section Sec)
433 {
434  const doublereal* pv1 = v1.pGetVec();
435  const doublereal* pv2 = v2.pGetVec();
436  const doublereal* pv3 = v3.pGetVec();
437 
438  return Vec3(
439  pv1[0]*dN3[Sec][0] + pv2[0]*dN3[Sec][1] + pv3[0]*dN3[Sec][2],
440  pv1[1]*dN3[Sec][0] + pv2[1]*dN3[Sec][1] + pv3[1]*dN3[Sec][2],
441  pv1[2]*dN3[Sec][0] + pv2[2]*dN3[Sec][1] + pv3[2]*dN3[Sec][2]);
442 }
443 
444 
445 Vec3
446 Beam::InterpDeriv(const Vec3& v1,
447  const Vec3& v2,
448  const Vec3& v3,
449  enum Section Sec)
450 {
451  const doublereal* pv1 = v1.pGetVec();
452  const doublereal* pv2 = v2.pGetVec();
453  const doublereal* pv3 = v3.pGetVec();
454 
455  return Vec3(
456  (pv1[0]*dN3P[Sec][0] + pv2[0]*dN3P[Sec][1]
457  + pv3[0]*dN3P[Sec][2])*dsdxi[Sec],
458  (pv1[1]*dN3P[Sec][0] + pv2[1]*dN3P[Sec][1]
459  + pv3[1]*dN3P[Sec][2])*dsdxi[Sec],
460  (pv1[2]*dN3P[Sec][0] + pv2[2]*dN3P[Sec][1]
461  + pv3[2]*dN3P[Sec][2])*dsdxi[Sec]);
462 }
463 
464 
465 void
466 Beam::DsDxi(void)
467 {
468  /* Validazione dati */
469  ASSERT(pNode[NODE1] != NULL);
470  ASSERT(pNode[NODE1]->GetNodeType() == Node::STRUCTURAL);
471  ASSERT(pNode[NODE2] != NULL);
472  ASSERT(pNode[NODE2]->GetNodeType() == Node::STRUCTURAL);
473  ASSERT(pNode[NODE3] != NULL);
474  ASSERT(pNode[NODE3]->GetNodeType() == Node::STRUCTURAL);
475 
476  /* Calcola il ds/dxi e le deformazioni iniziali */
477  Mat3x3 RNod[NUMNODES];
478  Vec3 xTmp[NUMNODES];
479  for (unsigned int i = 0; i < NUMNODES; i++) {
480  RNod[i] = pNode[i]->GetRCurr();
481  xTmp[i] = pNode[i]->GetXCurr() + RNod[i]*f[i];
482  }
483 
484  Vec3 xGrad[NUMSEZ];
485  for (unsigned int i = 0; i < NUMSEZ; i++) {
486  /* temporary */
487  dsdxi[i] = 1.;
488  xGrad[i] = InterpDeriv(xTmp[NODE1],
489  xTmp[NODE2],
490  xTmp[NODE3],
491  Beam::Section(i));
492  doublereal d = xGrad[i].Dot();
493  if (d > std::numeric_limits<doublereal>::epsilon()) {
494  /* final */
495  dsdxi[i] = 1./std::sqrt(d);
496 
497  } else {
498  silent_cerr("warning, beam " << GetLabel()
499  << " has singular metric; aborting..." << std::endl);
500 
501  throw ErrNullNorm(MBDYN_EXCEPT_ARGS);
502  }
503  }
504 
505  /* Calcola le deformazioni iniziali */
506  for (unsigned int i = 0; i < NUMSEZ; i++) {
507  L0[i] = R[i].MulTV(InterpDeriv(xTmp[NODE1],
508  xTmp[NODE2],
509  xTmp[NODE3],
510  Beam::Section(i)));
511  pD[i]->Update(Zero6);
512  DRef[i] = MultRMRt(pD[i]->GetFDE(), R[i]);
513  }
514 }
515 
516 
517 /* Calcola la velocita' angolare delle sezioni a partire da quelle dei nodi */
518 void
519 Beam::Omega0(void)
520 {
521  /* Validazione dati */
522  ASSERT(pNode[NODE1] != NULL);
523  ASSERT(pNode[NODE1]->GetNodeType() == Node::STRUCTURAL);
524  ASSERT(pNode[NODE2] != NULL);
525  ASSERT(pNode[NODE2]->GetNodeType() == Node::STRUCTURAL);
526  ASSERT(pNode[NODE3] != NULL);
527  ASSERT(pNode[NODE3]->GetNodeType() == Node::STRUCTURAL);
528 
529  /* Modo consistente: */
530  Mat3x3 RNod[NUMNODES];
531  Vec3 w[NUMNODES];
532  for (unsigned int i = 0; i < NUMNODES; i++) {
533  RNod[i] = pNode[i]->GetRCurr()*RNode[i];
534  w[i] = pNode[i]->GetWCurr();
535  }
536 
537  /* Calcolo i parametri di rotazione dei nodi di estremo rispetto a quello
538  * centrale, nell'ipotesi (realistica) che siano limitati */
539  Vec3 g1(CGR_Rot::Param, RNod[NODE2].MulTM(RNod[NODE1]));
540  Vec3 g3(CGR_Rot::Param, RNod[NODE2].MulTM(RNod[NODE3]));
541 
542  /* Calcolo le derivate dei parametri di rotazione ai nodi di estremo; in
543  * quello centrale si assume che sia uguale alla velocita' di rotazione */
544  Vec3 g1P(Mat3x3(CGR_Rot::MatGm1, g1)*w[NODE1]);
545  Vec3 g3P(Mat3x3(CGR_Rot::MatGm1, g3)*w[NODE3]);
546 
547  for (unsigned int i = 0; i < NUMSEZ; i++) {
548  Vec3 gTmp(g1*dN3[i][NODE1]+g3*dN3[i][NODE3]);
549  Vec3 gPTmp(g1P*dN3[i][NODE1]+w[NODE2]*dN3[i][NODE2]+g3P*dN3[i][NODE3]);
550  Omega[i] = Mat3x3(CGR_Rot::MatG, gTmp)*gPTmp;
551  }
552 
553 #if 0
554  /* Modo brutale: interpolo le velocita' dei nodi */
555  Vec3 w[NUMNODES];
556  for (unsigned int i = 0; i < NUMNODES; i++) {
557  w[i] = pNode[i]->GetWCurr();
558  }
559  for (unsigned int i = 0; i < NUMSEZ; i++) {
560  Omega[i] = (w{NODE1]*dN3[i][NODE1]
561  + w[NODE2]*dN3[i][NODE2]
562  + w[NODE3]*dN3[i][NODE3]);
563  }
564 #endif /* 0 */
565 }
566 
567 
568 /* Contributo al file di restart */
569 std::ostream&
570 Beam::Restart(std::ostream& out) const
571 {
572  return Restart_(out)<< ';' << std::endl;
573 }
574 
575 std::ostream&
576 Beam::Restart_(std::ostream& out) const
577 {
578  out << " beam: " << GetLabel();
579  for (unsigned int i = 0; i < NUMNODES; i++) {
580  out << ", " << pNode[i]->GetLabel() << ", reference, node, ",
581  f[i].Write(out, ", ");
582  }
583 
584  for (unsigned int i = 0; i < NUMSEZ; i++) {
585  out << ", reference, global, 1, ",
586  (R[i].GetVec(1)).Write(out, ", ")
587  << ", 2, ", (R[i].GetVec(2)).Write(out, ", ") << ", ",
588  pD[i]->pGetConstLaw()->Restart(out);
589  }
590 
591  return out;
592 }
593 
594 void
595 Beam::AfterConvergence(const VectorHandler& X, const VectorHandler& XP)
596 {
597  for (unsigned i = 0; i < NUMSEZ; i++) {
598  RPrev[i] = R[i];
599  DefLocPrev[i] = DefLoc[i];
600  pD[i]->AfterConvergence(DefLoc[i]);
601  }
602 }
603 
604 /* Inverse Dynamics: */
605 void
606 Beam::AfterConvergence(const VectorHandler& X, const VectorHandler& XP, const VectorHandler& XPP)
607 {
608  AfterConvergence(X, XP);
609 }
610 
611 /* Assembla la matrice */
612 void
613 Beam::AssStiffnessMat(FullSubMatrixHandler& WMA,
614  FullSubMatrixHandler& /* WMB */ ,
615  doublereal dCoef,
616  const VectorHandler& /* XCurr */ ,
617  const VectorHandler& /* XPrimeCurr */ )
618 {
619  DEBUGCOUTFNAME("Beam::AssStiffnessMat");
620 
621  /* La matrice arriva gia' dimensionata e con gli indici di righe e colonne
622  * a posto */
623 
624  /* offset nel riferimento globale */
625  Vec3 fTmp[NUMNODES];
626  for (unsigned int i = 0; i < NUMNODES; i++) {
627  fTmp[i] = pNode[i]->GetRCurr()*f[i];
628  }
629 
630  Mat6x6 AzTmp[NUMSEZ][NUMNODES];
631 
632  /* per ogni punto di valutazione: */
633  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
634  for (unsigned int i = 0; i < NUMNODES; i++) {
635  /* Delta - deformazioni */
636  AzTmp[iSez][i] = Mat6x6(mb_deye<Mat3x3>(dN3P[iSez][i]*dsdxi[iSez]*dCoef),
637  Zero3x3,
638  Mat3x3(MatCross, L[iSez]*(dN3[iSez][i]*dCoef) - fTmp[i]*(dN3P[iSez][i]*dsdxi[iSez]*dCoef)),
639  mb_deye<Mat3x3>(dN3P[iSez][i]*dsdxi[iSez]*dCoef));
640  /* Delta - azioni interne */
641  AzTmp[iSez][i] = DRef[iSez]*AzTmp[iSez][i];
642  /* Correggo per la rotazione da locale a globale */
643  AzTmp[iSez][i].SubMat12(Mat3x3(MatCross, Az[iSez].GetVec1()*(dN3[iSez][i]*dCoef)));
644  AzTmp[iSez][i].SubMat22(Mat3x3(MatCross, Az[iSez].GetVec2()*(dN3[iSez][i]*dCoef)));
645  }
646  } /* end ciclo sui punti di valutazione */
647 
648  Vec3 bTmp[2];
649 
650  /* ciclo sulle equazioni */
651  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
652  bTmp[0] = p[iSez] - pNode[iSez]->GetXCurr();
653  bTmp[1] = p[iSez] - pNode[iSez + 1]->GetXCurr();
654 
655  unsigned int iRow1 = iSez*6 + 1;
656  unsigned int iRow2 = iRow1 + 6;
657 
658  for (unsigned int i = 0; i < NUMNODES; i++) {
659  /* Equazione all'indietro: */
660  WMA.Sub(iRow1, 6*i + 1, AzTmp[iSez][i].GetMat11());
661  WMA.Sub(iRow1, 6*i + 4, AzTmp[iSez][i].GetMat12());
662 
663  WMA.Sub(iRow1 + 3, 6*i + 1,
664  AzTmp[iSez][i].GetMat21()
665  - Mat3x3(MatCross, Az[iSez].GetVec1()*(dCoef*dN3[iSez][i]))
666  + bTmp[0].Cross(AzTmp[iSez][i].GetMat11()));
667  WMA.Sub(iRow1 + 3, 6*i + 4,
668  AzTmp[iSez][i].GetMat22()
669  - Mat3x3(MatCrossCross, Az[iSez].GetVec1()*(-dCoef*dN3[iSez][i]), fTmp[i])
670  + bTmp[0].Cross(AzTmp[iSez][i].GetMat12()));
671 
672  /* Equazione in avanti: */
673  WMA.Add(iRow2, 6*i+1, AzTmp[iSez][i].GetMat11());
674  WMA.Add(iRow2, 6*i+4, AzTmp[iSez][i].GetMat12());
675 
676  WMA.Add(iRow2 + 3, 6*i + 1,
677  AzTmp[iSez][i].GetMat21()
678  - Mat3x3(MatCross, Az[iSez].GetVec1()*(dCoef*dN3[iSez][i]))
679  + bTmp[1].Cross(AzTmp[iSez][i].GetMat11()));
680  WMA.Add(iRow2 + 3, 6*i + 4,
681  AzTmp[iSez][i].GetMat22()
682  + Mat3x3(MatCrossCross, Az[iSez].GetVec1()*(dCoef*dN3[iSez][i]), fTmp[i])
683  + bTmp[1].Cross(AzTmp[iSez][i].GetMat12()));
684  }
685 
686  /* correzione alle equazioni */
687  Mat3x3 FTmp(MatCross, Az[iSez].GetVec1()*dCoef);
688  WMA.Sub(iRow1 + 3, 6*iSez + 1, FTmp);
689  WMA.Add(iRow2 + 3, 6*iSez + 7, FTmp);
690 
691  } /* end ciclo sui punti di valutazione */
692 }
693 
694 
695 /* Assembla il residuo */
696 void
697 Beam::AssStiffnessVec(SubVectorHandler& WorkVec,
698  doublereal /* dCoef */ ,
699  const VectorHandler& /* XCurr */ ,
700  const VectorHandler& /* XPrimeCurr */ )
701 {
702  DEBUGCOUTFNAME("Beam::AssStiffnessVec");
703 
704  /* Riceve il vettore gia' dimensionato e con gli indici a posto
705  * per scrivere il residuo delle equazioni di equilibrio dei tre nodi */
706 
707  /* Per la trattazione teorica, il riferimento e' il file ul-travi.tex
708  * (ora e' superato) */
709 
710  if (bFirstRes) {
711  bFirstRes = false; /* AfterPredict ha gia' calcolato tutto */
712  } else {
713  Vec3 gNod[NUMNODES];
714  Vec3 xTmp[NUMNODES];
715 
716  for (unsigned int i = 0; i < NUMNODES; i++) {
717  gNod[i] = pNode[i]->GetgCurr();
718  xTmp[i] = pNode[i]->GetXCurr() + pNode[i]->GetRCurr()*f[i];
719  }
720 
721  Mat3x3 RDelta[NUMSEZ];
722  Vec3 gGrad[NUMSEZ];
723 
724  /* Aggiorna le grandezze della trave nei punti di valutazione */
725  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
726 
727  /* Posizione */
728  p[iSez] = InterpState(xTmp[NODE1], xTmp[NODE2], xTmp[NODE3], Beam::Section(iSez));
729 
730  /* Matrici di rotazione */
731  g[iSez] = InterpState(gNod[NODE1], gNod[NODE2], gNod[NODE3], Beam::Section(iSez));
732  RDelta[iSez] = Mat3x3(CGR_Rot::MatR, g[iSez]);
733  R[iSez] = RDelta[iSez]*RRef[iSez];
734 
735  /* Derivate della posizione */
736  L[iSez] = InterpDeriv(xTmp[NODE1], xTmp[NODE2], xTmp[NODE3], Beam::Section(iSez));
737 
738  /* Derivate dei parametri di rotazione */
739  gGrad[iSez] = InterpDeriv(gNod[NODE1], gNod[NODE2], gNod[NODE3], Beam::Section(iSez));
740 
741  /* Calcola le deformazioni nel sistema locale nei punti di valutazione */
742  DefLoc[iSez] = Vec6(R[iSez].MulTV(L[iSez]) - L0[iSez],
743  R[iSez].MulTV(Mat3x3(CGR_Rot::MatG, g[iSez])*gGrad[iSez]) + DefLocRef[iSez].GetVec2());
744 
745  /* Calcola le azioni interne */
746  pD[iSez]->Update(DefLoc[iSez]);
747  AzLoc[iSez] = pD[iSez]->GetF();
748 
749  /* corregge le azioni interne locali (piezo, ecc) */
750  AddInternalForces(AzLoc[iSez], iSez);
751 
752  /* Porta le azioni interne nel sistema globale */
753  Az[iSez] = MultRV(AzLoc[iSez], R[iSez]);
754  }
755  }
756 
757  WorkVec.Add(1, Az[S_I].GetVec1());
758  WorkVec.Add(4, (p[S_I] - pNode[NODE1]->GetXCurr()).Cross(Az[S_I].GetVec1())
759  + Az[S_I].GetVec2());
760  WorkVec.Add(7, Az[SII].GetVec1() - Az[S_I].GetVec1());
761  WorkVec.Add(10, Az[SII].GetVec2() - Az[S_I].GetVec2()
762  + (p[SII] - pNode[NODE2]->GetXCurr()).Cross(Az[SII].GetVec1())
763  - (p[S_I] - pNode[NODE2]->GetXCurr()).Cross(Az[S_I].GetVec1()));
764  WorkVec.Sub(13, Az[SII].GetVec1());
765  WorkVec.Add(16, Az[SII].GetVec1().Cross(p[SII] - pNode[NODE3]->GetXCurr())
766  - Az[SII].GetVec2());
767 }
768 
769 
770 /* assemblaggio jacobiano */
771 VariableSubMatrixHandler&
772 Beam::AssJac(VariableSubMatrixHandler& WorkMat,
773  doublereal dCoef,
774  const VectorHandler& XCurr,
775  const VectorHandler& XPrimeCurr)
776 {
777  DEBUGCOUTFNAME("Beam::AssJac => AssStiffnessMat");
778 
779  integer iNode1FirstMomIndex = pNode[NODE1]->iGetFirstMomentumIndex();
780  integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
781  integer iNode2FirstMomIndex = pNode[NODE2]->iGetFirstMomentumIndex();
782  integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
783  integer iNode3FirstMomIndex = pNode[NODE3]->iGetFirstMomentumIndex();
784  integer iNode3FirstPosIndex = pNode[NODE3]->iGetFirstPositionIndex();
785 
786  FullSubMatrixHandler& WM = WorkMat.SetFull();
787 
788  /* Dimensiona la matrice, la azzera e pone gli indici corretti */
789  if (bConsistentInertia) {
790  WM.ResizeReset(36, 18);
791  } else {
792  WM.ResizeReset(18, 18);
793  }
794 
795  for (int iCnt = 1; iCnt <= 6; iCnt++) {
796  WM.PutRowIndex(iCnt, iNode1FirstMomIndex + iCnt);
797  WM.PutColIndex(iCnt, iNode1FirstPosIndex + iCnt);
798  WM.PutRowIndex(6 + iCnt, iNode2FirstMomIndex + iCnt);
799  WM.PutColIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
800  WM.PutRowIndex(12 + iCnt, iNode3FirstMomIndex + iCnt);
801  WM.PutColIndex(12 + iCnt, iNode3FirstPosIndex + iCnt);
802  }
803 
804  AssStiffnessMat(WM, WM, dCoef, XCurr, XPrimeCurr);
805 
806  if (bConsistentInertia) {
807  for (int iCnt = 1; iCnt <= 6; iCnt++) {
808  WM.PutRowIndex(18 + iCnt, iNode1FirstPosIndex + iCnt);
809  WM.PutRowIndex(24 + iCnt, iNode2FirstPosIndex + iCnt);
810  WM.PutRowIndex(30 + iCnt, iNode3FirstPosIndex + iCnt);
811  }
812 
813  AssInertiaMat(WM, WM, dCoef, XCurr, XPrimeCurr);
814  }
815 
816  return WorkMat;
817 }
818 
819 
820 /* assemblaggio matrici per autovalori */
821 void
822 Beam::AssMats(VariableSubMatrixHandler& WorkMatA,
823  VariableSubMatrixHandler& WorkMatB,
824  const VectorHandler& XCurr,
825  const VectorHandler& XPrimeCurr)
826 {
827  DEBUGCOUTFNAME("Beam::AssMats => AssStiffnessMat");
828 
829  integer iNode1FirstMomIndex = pNode[NODE1]->iGetFirstMomentumIndex();
830  integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
831  integer iNode2FirstMomIndex = pNode[NODE2]->iGetFirstMomentumIndex();
832  integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
833  integer iNode3FirstMomIndex = pNode[NODE3]->iGetFirstMomentumIndex();
834  integer iNode3FirstPosIndex = pNode[NODE3]->iGetFirstPositionIndex();
835 
836  FullSubMatrixHandler& WMA = WorkMatA.SetFull();
837  FullSubMatrixHandler& WMB = WorkMatB.SetFull();
838 
839  /* Dimensiona la matrice, la azzera e pone gli indici corretti */
840  if (bConsistentInertia) {
841  WMA.ResizeReset(36, 18);
842  WMB.ResizeReset(36, 18);
843  } else {
844  WMA.ResizeReset(18, 18);
845  WorkMatB.SetNullMatrix();
846  }
847 
848  for (int iCnt = 1; iCnt <= 6; iCnt++) {
849  WMA.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
850  WMA.PutColIndex(iCnt, iNode1FirstPosIndex+iCnt);
851  WMA.PutRowIndex(6+iCnt, iNode2FirstMomIndex+iCnt);
852  WMA.PutColIndex(6+iCnt, iNode2FirstPosIndex+iCnt);
853  WMA.PutRowIndex(12+iCnt, iNode3FirstMomIndex+iCnt);
854  WMA.PutColIndex(12+iCnt, iNode3FirstPosIndex+iCnt);
855  }
856 
857  AssStiffnessMat(WMA, WMA, 1., XCurr, XPrimeCurr);
858 
859  if (bConsistentInertia) {
860  for (unsigned int iCnt = 1; iCnt <= 6; iCnt++) {
861  WMB.PutRowIndex(iCnt, iNode1FirstMomIndex + iCnt);
862  WMB.PutColIndex(iCnt, iNode1FirstPosIndex + iCnt);
863  WMB.PutRowIndex(6 + iCnt, iNode2FirstMomIndex + iCnt);
864  WMB.PutColIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
865  WMB.PutRowIndex(12 + iCnt, iNode3FirstMomIndex + iCnt);
866  WMB.PutColIndex(12 + iCnt, iNode3FirstPosIndex + iCnt);
867  }
868 
869  for (unsigned int iCnt = 1; iCnt <= 6; iCnt++) {
870  WMA.PutRowIndex(18 + iCnt, iNode1FirstPosIndex + iCnt);
871  WMA.PutRowIndex(24 + iCnt, iNode2FirstPosIndex + iCnt);
872  WMA.PutRowIndex(30 + iCnt, iNode3FirstPosIndex + iCnt);
873 
874  WMB.PutRowIndex(18 + iCnt, iNode1FirstPosIndex + iCnt);
875  WMB.PutRowIndex(24 + iCnt, iNode2FirstPosIndex + iCnt);
876  WMB.PutRowIndex(30 + iCnt, iNode3FirstPosIndex + iCnt);
877  }
878 
879  AssInertiaMat(WMA, WMB, 1., XCurr, XPrimeCurr);
880  }
881 }
882 
883 
884 /* assemblaggio residuo */
885 SubVectorHandler&
886 Beam::AssRes(SubVectorHandler& WorkVec,
887  doublereal dCoef,
888  const VectorHandler& XCurr,
889  const VectorHandler& XPrimeCurr)
890 {
891  DEBUGCOUTFNAME("Beam::AssRes => AssStiffnessVec");
892 
893  integer iNode1FirstMomIndex = pNode[NODE1]->iGetFirstMomentumIndex();
894  integer iNode2FirstMomIndex = pNode[NODE2]->iGetFirstMomentumIndex();
895  integer iNode3FirstMomIndex = pNode[NODE3]->iGetFirstMomentumIndex();
896 
897  /* Dimensiona il vettore, lo azzera e pone gli indici corretti */
898  if (bConsistentInertia) {
899  WorkVec.ResizeReset(36);
900  } else {
901  WorkVec.ResizeReset(18);
902  }
903 
904  for (unsigned int iCnt = 1; iCnt <= 6; iCnt++) {
905  WorkVec.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
906  WorkVec.PutRowIndex(6 + iCnt, iNode2FirstMomIndex + iCnt);
907  WorkVec.PutRowIndex(12 + iCnt, iNode3FirstMomIndex + iCnt);
908  }
909 
910  AssStiffnessVec(WorkVec, dCoef, XCurr, XPrimeCurr);
911 
912  if (bConsistentInertia) {
913  integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
914  integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
915  integer iNode3FirstPosIndex = pNode[NODE3]->iGetFirstPositionIndex();
916 
917  for (unsigned int iCnt = 1; iCnt <= 6; iCnt++) {
918  WorkVec.PutRowIndex(18 + iCnt, iNode1FirstPosIndex + iCnt);
919  WorkVec.PutRowIndex(24 + iCnt, iNode2FirstPosIndex + iCnt);
920  WorkVec.PutRowIndex(30 + iCnt, iNode3FirstPosIndex + iCnt);
921  }
922 
923  AssInertiaVec(WorkVec, dCoef, XCurr, XPrimeCurr);
924  }
925 
926  return WorkVec;
927 }
928 
929 /* Inverse Dynamics residual assembly */
930 SubVectorHandler&
931 Beam::AssRes(SubVectorHandler& WorkVec,
932  const VectorHandler& XCurr ,
933  const VectorHandler& XPrimeCurr ,
934  const VectorHandler& XPrimePrimeCurr ,
935  InverseDynamics::Order iOrder)
936 {
937  DEBUGCOUTFNAME("Beam::AssRes => AssStiffnessVec");
938 
939  integer iNode1FirstMomIndex = pNode[NODE1]->iGetFirstPositionIndex();
940  integer iNode2FirstMomIndex = pNode[NODE2]->iGetFirstPositionIndex();
941  integer iNode3FirstMomIndex = pNode[NODE3]->iGetFirstPositionIndex();
942 
943  /* Dimensiona il vettore, lo azzera e pone gli indici corretti */
944  WorkVec.ResizeReset(18);
945 
946  for (unsigned int iCnt = 1; iCnt <= 6; iCnt++) {
947  WorkVec.PutRowIndex(iCnt, iNode1FirstMomIndex+iCnt);
948  WorkVec.PutRowIndex(6 + iCnt, iNode2FirstMomIndex + iCnt);
949  WorkVec.PutRowIndex(12 + iCnt, iNode3FirstMomIndex + iCnt);
950  }
951 
952  AssStiffnessVec(WorkVec, 1., XCurr, XPrimeCurr);
953 
954  return WorkVec;
955 }
956 
957 /* inverse dynamics capable element */
958 bool
959 Beam::bInverseDynamics(void) const
960 {
961  return true;
962 }
963 
964 /* Settings iniziali, prima della prima soluzione */
965 void
966 Beam::SetValue(DataManager *pDM,
967  VectorHandler& /* X */ , VectorHandler& /* XP */ ,
968  SimulationEntity::Hints *ph)
969 {
970  /* Aggiorna le grandezze della trave nei punti di valutazione */
971  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
972  RRef[iSez] = R[iSez];
973  LRef[iSez] = L[iSez];
974  DefLocRef[iSez] = DefLoc[iSez];
975  AzRef[iSez] = Az[iSez];
976 
977  /* Aggiorna il legame costitutivo di riferimento
978  * (la deformazione e' gia' stata aggiornata dall'ultimo residuo) */
979  DRef[iSez] = MultRMRt(pD[iSez]->GetFDE(), RRef[iSez]);
980  }
981 
982  bFirstRes = true;
983 }
984 
985 
986 /* Prepara i parametri di riferimento dopo la predizione */
987 void
988 Beam::AfterPredict(VectorHandler& /* X */ ,
989  VectorHandler& /* XP */ )
990 {
991  /* Calcola le deformazioni, aggiorna il legame costitutivo e crea la FDE */
992 
993  /* Recupera i dati dei nodi */
994  Vec3 gNod[NUMNODES];
995  Vec3 xTmp[NUMNODES];
996 
997  for (unsigned int i = 0; i < NUMNODES; i++) {
998  gNod[i] = pNode[i]->GetgRef();
999  xTmp[i] = pNode[i]->GetXCurr() + pNode[i]->GetRRef()*f[i];
1000  }
1001 
1002  Mat3x3 RDelta[NUMSEZ];
1003  Vec3 gGrad[NUMSEZ];
1004 
1005  /* Aggiorna le grandezze della trave nei punti di valutazione */
1006  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1007 
1008  /* Posizione */
1009  p[iSez] = InterpState(xTmp[NODE1], xTmp[NODE2], xTmp[NODE3], Beam::Section(iSez));
1010 
1011  /* Matrici di rotazione */
1012  g[iSez] = InterpState(gNod[NODE1], gNod[NODE2], gNod[NODE3], Beam::Section(iSez));
1013  RDelta[iSez] = Mat3x3(CGR_Rot::MatR, g[iSez]);
1014  R[iSez] = RRef[iSez] = RDelta[iSez]*RPrev[iSez];
1015 
1016  /* Derivate della posizione */
1017  L[iSez] = LRef[iSez]
1018  = InterpDeriv(xTmp[NODE1], xTmp[NODE2], xTmp[NODE3], Beam::Section(iSez));
1019 
1020  /* Derivate dei parametri di rotazione */
1021  gGrad[iSez]
1022  = InterpDeriv(gNod[NODE1], gNod[NODE2], gNod[NODE3], Beam::Section(iSez));
1023 
1024  /* Calcola le deformazioni nel sistema locale nei punti di valutazione */
1025  DefLoc[iSez] = DefLocRef[iSez]
1026  = Vec6(R[iSez].MulTV(L[iSez]) - L0[iSez],
1027  R[iSez].MulTV(Mat3x3(CGR_Rot::MatG, g[iSez])*gGrad[iSez]) + DefLocPrev[iSez].GetVec2());
1028 
1029  /* Calcola le azioni interne */
1030  pD[iSez]->Update(DefLoc[iSez]);
1031  AzLoc[iSez] = pD[iSez]->GetF();
1032 
1033  /* corregge le azioni interne locali (piezo, ecc) */
1034  AddInternalForces(AzLoc[iSez], iSez);
1035 
1036  /* Porta le azioni interne nel sistema globale */
1037  Az[iSez] = AzRef[iSez] = MultRV(AzLoc[iSez], R[iSez]);
1038 
1039  /* Aggiorna il legame costitutivo di riferimento */
1040  DRef[iSez] = MultRMRt(pD[iSez]->GetFDE(), RRef[iSez]);
1041  }
1042 
1043  bFirstRes = true;
1044 }
1045 
1046 void
1047 Beam::OutputPrepare(OutputHandler &OH)
1048 {
1049  if (bToBeOutput()) {
1050 #ifdef USE_NETCDF
1051  if (OH.UseNetCDF(OutputHandler::BEAMS)) {
1052  ASSERT(OH.IsOpen(OutputHandler::NETCDF));
1053 
1054  const char *type = 0;
1055  switch (GetBeamType()) {
1056  case Beam::ELASTIC:
1057  type = "elastic";
1058  break;
1059 
1060  case Beam::VISCOELASTIC:
1061  type = "viscoelastic";
1062  break;
1063 
1064  case Beam::PIEZOELECTRICELASTIC:
1065  type = "piezoelectric elastic";
1066  break;
1067 
1068  case Beam::PIEZOELECTRICVISCOELASTIC:
1069  type = "piezoelectric viscoelastic";
1070  break;
1071 
1072  default:
1073  type = "unknown";
1074  break;
1075  }
1076 
1077  std::ostringstream os;
1078  os << "elem.beam." << GetLabel();
1079 
1080  (void)OH.CreateVar(os.str(), type);
1081 
1082  os << '.';
1083  std::string name(os.str());
1084 
1085  unsigned uOutputFlags = (fToBeOutput() & ToBeOutput::OUTPUT_PRIVATE_MASK);
1086 
1087  static const char *sez[] = { "I", "II" };
1088  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1089  os.str("");
1090  os << "evaluation point " << sez[iSez] << " ";
1091  std::string ep(os.str());
1092 
1093  if (uOutputFlags & Beam::OUTPUT_EP_X) {
1094  Var_X[iSez] = OH.CreateVar<Vec3>(name + "X_" + sez[iSez], "m",
1095  ep + "global position vector (X, Y, Z)");
1096  }
1097 
1098  if (uOutputFlags & Beam::OUTPUT_EP_R) {
1099  Var_Phi[iSez] = OH.CreateRotationVar(name,
1100  std::string("_") + sez[iSez], od, ep + "global");
1101  }
1102 
1103  if (uOutputFlags & Beam::OUTPUT_EP_F) {
1104  Var_F[iSez] = OH.CreateVar<Vec3>(name + "F_" + sez[iSez], "N",
1105  ep + "internal force in local frame (F_X, F_Y, F_Z)");
1106  }
1107 
1108  if (uOutputFlags & Beam::OUTPUT_EP_M) {
1109  Var_M[iSez] = OH.CreateVar<Vec3>(name + "M_" + sez[iSez], "Nm",
1110  ep + "internal moment in local frame (M_X, M_Y, M_Z)");
1111  }
1112 
1113  if (uOutputFlags & Beam::OUTPUT_EP_NU) {
1114  Var_Nu[iSez] = OH.CreateVar<Vec3>(name + "nu_" + sez[iSez], "-",
1115  ep + "linear strain in local frame (nu_X, nu_Y, nu_Z)");
1116  }
1117 
1118  if (uOutputFlags & Beam::OUTPUT_EP_K) {
1119  Var_K[iSez] = OH.CreateVar<Vec3>(name + "k_" + sez[iSez], "1/m",
1120  ep + "angular strain in local frame (K_X, K_Y, K_Z)");
1121  }
1122 
1123  if (uOutputFlags & Beam::OUTPUT_EP_NUP) {
1124  Var_NuP[iSez] = OH.CreateVar<Vec3>(name + "nuP_" + sez[iSez], "1/s",
1125  ep + "linear strain rate in local frame (nuP_X, nuP_Y, nuP_Z)");
1126  }
1127 
1128  if (uOutputFlags & Beam::OUTPUT_EP_KP) {
1129  Var_KP[iSez] = OH.CreateVar<Vec3>(name + "kP_" + sez[iSez], "1/ms",
1130  ep + "angular strain rate in local frame (KP_X, KP_Y, KP_Z)");
1131  }
1132  }
1133  }
1134 #endif // USE_NETCDF
1135  }
1136 }
1137 
1138 /* output; si assume che ogni tipo di elemento sappia, attraverso
1139  * l'OutputHandler, dove scrivere il proprio output */
1140 void
1141 Beam::Output(OutputHandler& OH) const
1142 {
1143  if (bToBeOutput()) {
1144 #ifdef USE_NETCDF
1145  if (OH.UseNetCDF(OutputHandler::BEAMS)) {
1146  for (unsigned iSez = 0; iSez < NUMSEZ; iSez++) {
1147  if (Var_X[iSez]) {
1148  Var_X[iSez]->put_rec(p[iSez].pGetVec(), OH.GetCurrentStep());
1149  }
1150 
1151  if (Var_Phi[iSez]) {
1152  Vec3 E;
1153  switch (od) {
1154  case EULER_123:
1155  E = MatR2EulerAngles123(R[iSez])*dRaDegr;
1156  break;
1157 
1158  case EULER_313:
1159  E = MatR2EulerAngles313(R[iSez])*dRaDegr;
1160  break;
1161 
1162  case EULER_321:
1163  E = MatR2EulerAngles321(R[iSez])*dRaDegr;
1164  break;
1165 
1166  case ORIENTATION_VECTOR:
1167  E = RotManip::VecRot(R[iSez]);
1168  break;
1169 
1170  case ORIENTATION_MATRIX:
1171  break;
1172 
1173  default:
1174  /* impossible */
1175  break;
1176  }
1177 
1178  switch (od) {
1179  case EULER_123:
1180  case EULER_313:
1181  case EULER_321:
1182  case ORIENTATION_VECTOR:
1183  Var_Phi[iSez]->put_rec(E.pGetVec(), OH.GetCurrentStep());
1184  break;
1185 
1186  case ORIENTATION_MATRIX:
1187  Var_Phi[iSez]->put_rec(R[iSez].pGetMat(), OH.GetCurrentStep());
1188  break;
1189 
1190  default:
1191  /* impossible */
1192  break;
1193  }
1194  }
1195 
1196  if (Var_F[iSez]) {
1197  Var_F[iSez]->put_rec(AzLoc[iSez].GetVec1().pGetVec(), OH.GetCurrentStep());
1198  }
1199 
1200  if (Var_M[iSez]) {
1201  Var_M[iSez]->put_rec(AzLoc[iSez].GetVec2().pGetVec(), OH.GetCurrentStep());
1202  }
1203 
1204  if (Var_Nu[iSez]) {
1205  Var_Nu[iSez]->put_rec(DefLoc[iSez].GetVec1().pGetVec(), OH.GetCurrentStep());
1206  }
1207 
1208  if (Var_K[iSez]) {
1209  Var_K[iSez]->put_rec(DefLoc[iSez].GetVec2().pGetVec(), OH.GetCurrentStep());
1210  }
1211 
1212  if (Var_NuP[iSez]) {
1213  Var_NuP[iSez]->put_rec(DefPrimeLoc[iSez].GetVec1().pGetVec(), OH.GetCurrentStep());
1214  }
1215 
1216  if (Var_KP[iSez]) {
1217  Var_KP[iSez]->put_rec(DefPrimeLoc[iSez].GetVec2().pGetVec(), OH.GetCurrentStep());
1218  }
1219  }
1220  }
1221 #endif /* USE_NETCDF */
1222 
1223  if (OH.UseText(OutputHandler::BEAMS)) {
1224  OH.Beams() << std::setw(8) << GetLabel()
1225  << " " << AzLoc[S_I].GetVec1()
1226  << " " << AzLoc[S_I].GetVec2()
1227  << " " << AzLoc[SII].GetVec1()
1228  << " " << AzLoc[SII].GetVec2()
1229  << std::endl;
1230  }
1231  }
1232 }
1233 
1234 /* Contributo allo jacobiano durante l'assemblaggio iniziale */
1235 VariableSubMatrixHandler&
1236 Beam::InitialAssJac(VariableSubMatrixHandler& WorkMat,
1237  const VectorHandler& XCurr)
1238 {
1239  DEBUGCOUTFNAME("Beam::InitialAssJac => AssStiffnessMat");
1240 
1241  /* Dimensiona la matrice, la azzera e pone gli indici corretti */
1242  FullSubMatrixHandler& WM = WorkMat.SetFull();
1243  WM.ResizeReset(18, 18);
1244 
1245  integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
1246  integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
1247  integer iNode3FirstPosIndex = pNode[NODE3]->iGetFirstPositionIndex();
1248 
1249  for (int iCnt = 1; iCnt <= 6; iCnt++) {
1250  WM.PutRowIndex(iCnt, iNode1FirstPosIndex + iCnt);
1251  WM.PutColIndex(iCnt, iNode1FirstPosIndex + iCnt);
1252  WM.PutRowIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
1253  WM.PutColIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
1254  WM.PutRowIndex(12 + iCnt, iNode3FirstPosIndex + iCnt);
1255  WM.PutColIndex(12 + iCnt, iNode3FirstPosIndex + iCnt);
1256  }
1257 
1258  AssStiffnessMat(WM, WM, 1., XCurr, XCurr);
1259 
1260  return WorkMat;
1261 }
1262 
1263 
1264 /* Contributo al residuo durante l'assemblaggio iniziale */
1265 SubVectorHandler& Beam::InitialAssRes(SubVectorHandler& WorkVec,
1266  const VectorHandler& XCurr)
1267 {
1268  DEBUGCOUTFNAME("Beam::InitialAssRes => AssStiffnessVec");
1269 
1270  /* Dimensiona il vettore, lo azzera e pone gli indici corretti */
1271  WorkVec.ResizeReset(18);
1272 
1273  integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
1274  integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
1275  integer iNode3FirstPosIndex = pNode[NODE3]->iGetFirstPositionIndex();
1276 
1277  for (int iCnt = 1; iCnt <= 6; iCnt++) {
1278  WorkVec.PutRowIndex(iCnt, iNode1FirstPosIndex + iCnt);
1279  WorkVec.PutRowIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
1280  WorkVec.PutRowIndex(12 + iCnt, iNode3FirstPosIndex + iCnt);
1281  }
1282 
1283  AssStiffnessVec(WorkVec, 1., XCurr, XCurr);
1284 
1285  return WorkVec;
1286 }
1287 
1288 
1289 const StructNode*
1290 Beam::pGetNode(unsigned int i) const
1291 {
1292  ASSERT(i >= 1 && i <= 3);
1293  switch (i) {
1294  case 1:
1295  case 2:
1296  case 3:
1297  return pNode[i-1];
1298  default:
1299  throw Beam::ErrGeneric(MBDYN_EXCEPT_ARGS);
1300  }
1301 }
1302 
1303 
1304 /* Beam - end */
1305 
1306 
1307 /* ViscoElasticBeam - begin */
1308 
1309 /* Costruttore normale */
1310 ViscoElasticBeam::ViscoElasticBeam(
1311  unsigned int uL,
1312  const StructNode* pN1,
1313  const StructNode* pN2,
1314  const StructNode* pN3,
1315  const Vec3& F1,
1316  const Vec3& F2,
1317  const Vec3& F3,
1318  const Mat3x3& R1,
1319  const Mat3x3& R2,
1320  const Mat3x3& R3,
1321  const Mat3x3& r_I,
1322  const Mat3x3& rII,
1323  const ConstitutiveLaw6D* pD_I,
1324  const ConstitutiveLaw6D* pDII,
1325  OrientationDescription ood,
1326  flag fOut)
1327 : Elem(uL, fOut),
1328 Beam(uL, pN1, pN2, pN3, F1, F2, F3, R1, R2, R3, r_I, rII, pD_I, pDII, ood, fOut)
1329 {
1330  Init();
1331 }
1332 
1333 
1334 /* Costruttore per la trave con forze d'inerzia consistenti */
1335 ViscoElasticBeam::ViscoElasticBeam(
1336  unsigned int uL,
1337  const StructNode* pN1,
1338  const StructNode* pN2,
1339  const StructNode* pN3,
1340  const Vec3& F1,
1341  const Vec3& F2,
1342  const Vec3& F3,
1343  const Mat3x3& R1,
1344  const Mat3x3& R2,
1345  const Mat3x3& R3,
1346  const Mat3x3& r_I, const Mat3x3& rII,
1347  const ConstitutiveLaw6D* pD_I,
1348  const ConstitutiveLaw6D* pDII,
1349  doublereal dM_I,
1350  const Vec3& s0_I, const Mat3x3& j0_I,
1351  doublereal dMII,
1352  const Vec3& s0II, const Mat3x3& j0II,
1353  OrientationDescription ood,
1354  flag fOut)
1355 : Elem(uL, fOut),
1356 Beam(uL, pN1, pN2, pN3, F1, F2, F3, R1, R2, R3, r_I, rII, pD_I, pDII,
1357  dM_I, s0_I, j0_I, dMII, s0II, j0II, ood, fOut)
1358 {
1359  Init();
1360 }
1361 
1362 void
1363 ViscoElasticBeam::Init(void)
1364 {
1365  gPrime[S_I] = gPrime[SII] = Zero3;
1366 
1367  LPrime[S_I] = LPrime[SII] = Zero3;
1368  LPrimeRef[S_I] = LPrimeRef[SII] = Zero3;
1369 
1370  DefPrimeLoc[S_I] = DefPrimeLoc[SII] = Zero6;
1371  DefPrimeLocRef[S_I] = DefPrimeLocRef[SII] = Zero6;
1372 
1373  /* Nota: DsDxi() viene chiamata dal costruttore di Beam */
1374  Beam::Omega0();
1375 
1376  OmegaRef[S_I] = Omega[S_I];
1377  OmegaRef[SII] = Omega[SII];
1378 
1379  const_cast<Mat6x6&>(ERef[S_I]) = MultRMRt(pD[S_I]->GetFDEPrime(), RRef[S_I]);
1380  const_cast<Mat6x6&>(ERef[SII]) = MultRMRt(pD[SII]->GetFDEPrime(), RRef[SII]);
1381 }
1382 
1383 void
1384 ViscoElasticBeam::AfterConvergence(const VectorHandler& X,
1385  const VectorHandler& XP)
1386 {
1387  for (unsigned i = 0; i < NUMSEZ; i++) {
1388  RPrev[i] = R[i];
1389  DefLocPrev[i] = DefLoc[i];
1390  pD[i]->AfterConvergence(DefLoc[i], DefPrimeLoc[i]);
1391  }
1392 }
1393 
1394 /* Inverse Dynamics */
1395 void
1396 ViscoElasticBeam::AfterConvergence(const VectorHandler& X,
1397  const VectorHandler& XP, const VectorHandler& XPP)
1398 {
1399  AfterConvergence(X, XP);
1400 }
1401 
1402 /* Assembla la matrice */
1403 void
1404 ViscoElasticBeam::AssStiffnessMat(FullSubMatrixHandler& WMA,
1405  FullSubMatrixHandler& WMB,
1406  doublereal dCoef,
1407  const VectorHandler& /* XCurr */ ,
1408  const VectorHandler& /* XPrimeCurr */ )
1409 {
1410  DEBUGCOUTFNAME("ViscoElasticBeam::AssStiffnessMat");
1411 
1412  /* La matrice arriva gia' dimensionata e con gli indici di righe e colonne
1413  * a posto */
1414 
1415  /* offset nel riferimento globale */
1416  Vec3 fTmp[NUMNODES];
1417  for (unsigned int i = 0; i < NUMNODES; i++) {
1418  fTmp[i] = pNode[i]->GetRCurr()*f[i];
1419  }
1420 
1421  Mat6x6 AzTmp[NUMSEZ][NUMNODES];
1422  Mat6x6 AzPrimeTmp[NUMSEZ][NUMNODES];
1423 
1424  /* per ogni punto di valutazione: */
1425  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1426  for (unsigned int i = 0; i < NUMNODES; i++) {
1427  /* Delta - deformazioni */
1428  AzTmp[iSez][i] = AzPrimeTmp[iSez][i]
1429  = Mat6x6(mb_deye<Mat3x3>(dN3P[iSez][i]*dsdxi[iSez]),
1430  Zero3x3,
1431  Mat3x3(MatCross, L[iSez]*dN3[iSez][i] - fTmp[i]*(dN3P[iSez][i]*dsdxi[iSez])),
1432  mb_deye<Mat3x3>(dN3P[iSez][i]*dsdxi[iSez]));
1433  AzTmp[iSez][i] = DRef[iSez]*AzTmp[iSez][i]*dCoef;
1434  AzTmp[iSez][i] +=
1435  ERef[iSez]*Mat6x6(Mat3x3(MatCross, Omega[iSez]*(-dN3P[iSez][i]*dsdxi[iSez]*dCoef)),
1436  Zero3x3,
1437  (Mat3x3(MatCross, LPrime[iSez]) - Mat3x3(MatCrossCross, Omega[iSez], L[iSez]))*(dN3[iSez][i]*dCoef)
1438  + Mat3x3(MatCrossCross, Omega[iSez], fTmp[i]*(dN3P[iSez][i]*dsdxi[iSez]*dCoef))
1439  + Mat3x3(MatCross, fTmp[i].Cross(pNode[i]->GetWCurr()*(dN3P[iSez][i]*dsdxi[iSez]*dCoef))),
1440  Mat3x3(MatCross, Omega[iSez]*(-dN3P[iSez][i]*dsdxi[iSez]*dCoef)));
1441  AzPrimeTmp[iSez][i] = ERef[iSez]*AzPrimeTmp[iSez][i];
1442 
1443  /* Correggo per la rotazione da locale a globale */
1444  AzTmp[iSez][i].SubMat12(Mat3x3(MatCross, Az[iSez].GetVec1()*(dN3[iSez][i]*dCoef)));
1445  AzTmp[iSez][i].SubMat22(Mat3x3(MatCross, Az[iSez].GetVec2()*(dN3[iSez][i]*dCoef)));
1446  }
1447  } /* end ciclo sui punti di valutazione */
1448 
1449 
1450  Vec3 bTmp[2];
1451 
1452  /* ciclo sulle equazioni */
1453  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1454  bTmp[0] = p[iSez] - pNode[iSez]->GetXCurr();
1455  bTmp[1] = p[iSez] - pNode[iSez + 1]->GetXCurr();
1456 
1457  unsigned int iRow1 = iSez*6 + 1;
1458  unsigned int iRow2 = iRow1 + 6;
1459 
1460  for (unsigned int i = 0; i < NUMNODES; i++) {
1461  /* Equazione all'indietro: */
1462  WMA.Sub(iRow1, 6*i + 1, AzTmp[iSez][i].GetMat11());
1463  WMA.Sub(iRow1, 6*i + 4, AzTmp[iSez][i].GetMat12());
1464 
1465  WMA.Sub(iRow1 + 3, 6*i + 1,
1466  AzTmp[iSez][i].GetMat21()
1467  - Mat3x3(MatCross, Az[iSez].GetVec1()*(dCoef*dN3[iSez][i]))
1468  + bTmp[0].Cross(AzTmp[iSez][i].GetMat11()));
1469  WMA.Sub(iRow1 + 3, 6*i + 4,
1470  AzTmp[iSez][i].GetMat22()
1471  - Mat3x3(MatCrossCross, Az[iSez].GetVec1()*(-dCoef*dN3[iSez][i]), fTmp[i])
1472  + bTmp[0].Cross(AzTmp[iSez][i].GetMat12()));
1473 
1474  /* Equazione in avanti: */
1475  WMA.Add(iRow2, 6*i + 1, AzTmp[iSez][i].GetMat11());
1476  WMA.Add(iRow2, 6*i + 4, AzTmp[iSez][i].GetMat12());
1477 
1478  WMA.Add(iRow2 + 3, 6*i + 1,
1479  AzTmp[iSez][i].GetMat21()
1480  - Mat3x3(MatCross, Az[iSez].GetVec1()*(dCoef*dN3[iSez][i]))
1481  + bTmp[1].Cross(AzTmp[iSez][i].GetMat11()));
1482  WMA.Add(iRow2 + 3, 6*i + 4,
1483  AzTmp[iSez][i].GetMat22()
1484  + Mat3x3(MatCrossCross, Az[iSez].GetVec1()*(dCoef*dN3[iSez][i]), fTmp[i])
1485  + bTmp[1].Cross(AzTmp[iSez][i].GetMat12()));
1486 
1487 
1488  /* Equazione viscosa all'indietro: */
1489  WMB.Sub(iRow1, 6*i + 1, AzPrimeTmp[iSez][i].GetMat11());
1490  WMB.Sub(iRow1, 6*i + 4, AzPrimeTmp[iSez][i].GetMat12());
1491 
1492  WMB.Sub(iRow1 + 3, 6*i + 1,
1493  AzPrimeTmp[iSez][i].GetMat21()
1494  + bTmp[0].Cross(AzPrimeTmp[iSez][i].GetMat11()));
1495  WMB.Sub(iRow1 + 3, 6*i + 4,
1496  AzPrimeTmp[iSez][i].GetMat22()
1497  + bTmp[0].Cross(AzPrimeTmp[iSez][i].GetMat12()));
1498 
1499  /* Equazione viscosa in avanti: */
1500  WMB.Add(iRow2, 6*i + 1, AzPrimeTmp[iSez][i].GetMat11());
1501  WMB.Add(iRow2, 6*i + 4, AzPrimeTmp[iSez][i].GetMat12());
1502 
1503  WMB.Add(iRow2 + 3, 6*i + 1,
1504  AzPrimeTmp[iSez][i].GetMat21()
1505  + bTmp[1].Cross(AzPrimeTmp[iSez][i].GetMat11()));
1506  WMB.Add(iRow2 + 3, 6*i + 4,
1507  AzPrimeTmp[iSez][i].GetMat22()
1508  + bTmp[1].Cross(AzPrimeTmp[iSez][i].GetMat12()));
1509  }
1510 
1511  /* correzione alle equazioni */
1512  Mat3x3 FTmp(MatCross, Az[iSez].GetVec1()*dCoef);
1513  WMA.Sub(iRow1 + 3, 6*iSez + 1, FTmp);
1514  WMA.Add(iRow2 + 3, 6*iSez + 7, FTmp);
1515 
1516  } /* end ciclo sui punti di valutazione */
1517 };
1518 
1519 
1520 /* Assembla il residuo */
1521 void
1522 ViscoElasticBeam::AssStiffnessVec(SubVectorHandler& WorkVec,
1523  doublereal /* dCoef */ ,
1524  const VectorHandler& /* XCurr */ ,
1525  const VectorHandler& /* XPrimeCurr */ )
1526 {
1527  DEBUGCOUTFNAME("ViscoElasticBeam::AssStiffnessVec");
1528 
1529  /* Riceve il vettore gia' dimensionato e con gli indici a posto
1530  * per scrivere il residuo delle equazioni di equilibrio dei tre nodi */
1531 
1532  /* Per la trattazione teorica, il riferimento e' il file ul-travi.tex
1533  * (ora e' superato) */
1534 
1535  if (bFirstRes) {
1536  bFirstRes = false; /* AfterPredict ha gia' calcolato tutto */
1537  } else {
1538  Vec3 gNod[NUMNODES];
1539  Vec3 xTmp[NUMNODES];
1540 
1541  Vec3 gPrimeNod[NUMNODES];
1542  Vec3 xPrimeTmp[NUMNODES];
1543 
1544  for (unsigned int i = 0; i < NUMNODES; i++) {
1545  gNod[i] = pNode[i]->GetgCurr();
1546  Vec3 fTmp = pNode[i]->GetRCurr()*f[i];
1547  xTmp[i] = pNode[i]->GetXCurr() + fTmp;
1548  gPrimeNod[i] = pNode[i]->GetgPCurr();
1549  xPrimeTmp[i] = pNode[i]->GetVCurr()+pNode[i]->GetWCurr().Cross(fTmp);
1550  }
1551 
1552  Mat3x3 RDelta[NUMSEZ];
1553  Vec3 gGrad[NUMSEZ];
1554  Vec3 gPrimeGrad[NUMSEZ];
1555 
1556  /* Aggiorna le grandezze della trave nei punti di valutazione */
1557  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1558 
1559  /* Posizione */
1560  p[iSez] = InterpState(xTmp[NODE1],
1561  xTmp[NODE2],
1562  xTmp[NODE3], Beam::Section(iSez));
1563 
1564  /* Matrici di rotazione */
1565  g[iSez] = InterpState(gNod[NODE1],
1566  gNod[NODE2],
1567  gNod[NODE3], Beam::Section(iSez));
1568  RDelta[iSez] = Mat3x3(CGR_Rot::MatR, g[iSez]);
1569  R[iSez] = RDelta[iSez]*RRef[iSez];
1570 
1571  /* Velocita' angolare della sezione */
1572  gPrime[iSez] = InterpState(gPrimeNod[NODE1],
1573  gPrimeNod[NODE2],
1574  gPrimeNod[NODE3], Beam::Section(iSez));
1575  Omega[iSez] = Mat3x3(CGR_Rot::MatG, g[iSez])*gPrime[iSez]
1576  + RDelta[iSez]*OmegaRef[iSez];
1577 
1578  /* Derivate della posizione */
1579  L[iSez] = InterpDeriv(xTmp[NODE1],
1580  xTmp[NODE2],
1581  xTmp[NODE3], Beam::Section(iSez));
1582 
1583  /* Derivate della velocita' */
1584  LPrime[iSez] = InterpDeriv(xPrimeTmp[NODE1],
1585  xPrimeTmp[NODE2],
1586  xPrimeTmp[NODE3], Beam::Section(iSez));
1587 
1588  /* Derivate dei parametri di rotazione */
1589  gGrad[iSez] = InterpDeriv(gNod[NODE1],
1590  gNod[NODE2],
1591  gNod[NODE3], Beam::Section(iSez));
1592 
1593  /* Derivate delle derivate spaziali dei parametri di rotazione */
1594  gPrimeGrad[iSez] = InterpDeriv(gPrimeNod[NODE1],
1595  gPrimeNod[NODE2],
1596  gPrimeNod[NODE3], Beam::Section(iSez));
1597 
1598  /* Calcola le deformazioni nel sistema locale nei punti di valutazione */
1599  DefLoc[iSez] = Vec6(R[iSez].MulTV(L[iSez]) - L0[iSez],
1600  R[iSez].MulTV(Mat3x3(CGR_Rot::MatG, g[iSez])*gGrad[iSez]) + DefLocRef[iSez].GetVec2());
1601 
1602  /* Calcola le velocita' di deformazione nel sistema locale nei punti di valutazione */
1603  DefPrimeLoc[iSez] = Vec6(R[iSez].MulTV(LPrime[iSez] + L[iSez].Cross(Omega[iSez])),
1604  R[iSez].MulTV(Mat3x3(CGR_Rot::MatG, g[iSez])*gPrimeGrad[iSez]
1605  + (Mat3x3(CGR_Rot::MatG, g[iSez])*gGrad[iSez]).Cross(Omega[iSez]))
1606  + DefPrimeLocRef[iSez].GetVec2());
1607 
1608  /* Calcola le azioni interne */
1609  pD[iSez]->Update(DefLoc[iSez], DefPrimeLoc[iSez]);
1610  AzLoc[iSez] = pD[iSez]->GetF();
1611 
1612  /* corregge le azioni interne locali (piezo, ecc) */
1613  AddInternalForces(AzLoc[iSez], iSez);
1614 
1615  /* Porta le azioni interne nel sistema globale */
1616  Az[iSez] = MultRV(AzLoc[iSez], R[iSez]);
1617  }
1618  }
1619 
1620  WorkVec.Add(1, Az[S_I].GetVec1());
1621  WorkVec.Add(4, (p[S_I] - pNode[NODE1]->GetXCurr()).Cross(Az[S_I].GetVec1())
1622  + Az[S_I].GetVec2());
1623  WorkVec.Add(7, Az[SII].GetVec1() - Az[S_I].GetVec1());
1624  WorkVec.Add(10, Az[SII].GetVec2() - Az[S_I].GetVec2()
1625  + (p[SII] - pNode[NODE2]->GetXCurr()).Cross(Az[SII].GetVec1())
1626  - (p[S_I] - pNode[NODE2]->GetXCurr()).Cross(Az[S_I].GetVec1()));
1627  WorkVec.Sub(13, Az[SII].GetVec1());
1628  WorkVec.Add(16, Az[SII].GetVec1().Cross(p[SII] - pNode[NODE3]->GetXCurr())
1629  - Az[SII].GetVec2());
1630 }
1631 
1632 
1633 /* Settings iniziali, prima della prima soluzione */
1634 void
1635 ViscoElasticBeam::SetValue(DataManager *pDM,
1636  VectorHandler& X, VectorHandler& XP,
1637  SimulationEntity::Hints *ph)
1638 {
1639  Beam::SetValue(pDM, X, XP, ph);
1640 
1641  /* Aggiorna le grandezze della trave nei punti di valutazione */
1642  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1643  OmegaRef[iSez] = Omega[iSez];
1644  LPrimeRef[iSez] = LPrime[iSez];
1645  DefPrimeLocRef[iSez] = DefPrimeLoc[iSez];
1646 
1647  /* Aggiorna il legame costitutivo di riferimento
1648  * (la deformazione e' gia' stata aggiornata dall'ultimo residuo) */
1649  ERef[iSez] = MultRMRt(pD[iSez]->GetFDEPrime(), RRef[iSez]);
1650  }
1651  ASSERT(bFirstRes == true);
1652 }
1653 
1654 
1655 /* Prepara i parametri di riferimento dopo la predizione */
1656 void
1657 ViscoElasticBeam::AfterPredict(VectorHandler& /* X */ ,
1658  VectorHandler& /* XP */ )
1659 {
1660  /* Calcola le deformazioni, aggiorna il legame costitutivo e crea la FDE */
1661 
1662  /* Recupera i dati dei nodi */
1663  Vec3 gNod[NUMNODES];
1664  Vec3 xTmp[NUMNODES];
1665 
1666  Vec3 gPrimeNod[NUMNODES];
1667  Vec3 xPrimeTmp[NUMNODES];
1668 
1669  for (unsigned int i = 0; i < NUMNODES; i++) {
1670  gNod[i] = pNode[i]->GetgRef();
1671  Vec3 fTmp = pNode[i]->GetRRef()*f[i];
1672  xTmp[i] = pNode[i]->GetXCurr() + fTmp;
1673  gPrimeNod[i] = pNode[i]->GetgPRef();
1674  xPrimeTmp[i] = pNode[i]->GetVCurr() + pNode[i]->GetWRef().Cross(fTmp);
1675  }
1676 
1677  Mat3x3 RDelta[NUMSEZ];
1678  Vec3 gGrad[NUMSEZ];
1679  Vec3 gPrimeGrad[NUMSEZ];
1680 
1681  /* Aggiorna le grandezze della trave nei punti di valutazione */
1682  for (unsigned int iSez = 0; iSez < NUMSEZ; iSez++) {
1683  /* Posizione */
1684  p[iSez] = InterpState(xTmp[NODE1],
1685  xTmp[NODE2],
1686  xTmp[NODE3], Beam::Section(iSez));
1687 
1688  /* Matrici di rotazione */
1689  g[iSez] = InterpState(gNod[NODE1],
1690  gNod[NODE2],
1691  gNod[NODE3], Beam::Section(iSez));
1692  RDelta[iSez] = Mat3x3(CGR_Rot::MatR, g[iSez]);
1693  R[iSez] = RRef[iSez] = RDelta[iSez]*RPrev[iSez];
1694 
1695  /* Velocita' angolare della sezione */
1696  gPrime[iSez] = InterpState(gPrimeNod[NODE1],
1697  gPrimeNod[NODE2],
1698  gPrimeNod[NODE3], Beam::Section(iSez));
1699  Omega[iSez] = OmegaRef[iSez] = Mat3x3(CGR_Rot::MatG, g[iSez])*gPrime[iSez];
1700 
1701  /* Derivate della posizione */
1702  L[iSez] = LRef[iSez] = InterpDeriv(xTmp[NODE1],
1703  xTmp[NODE2],
1704  xTmp[NODE3], Beam::Section(iSez));
1705 
1706  /* Derivate della velocita' */
1707  LPrime[iSez] = LPrimeRef[iSez]
1708  = InterpDeriv(xPrimeTmp[NODE1],
1709  xPrimeTmp[NODE2],
1710  xPrimeTmp[NODE3], Beam::Section(iSez));
1711 
1712  /* Derivate dei parametri di rotazione */
1713  gGrad[iSez] = InterpDeriv(gNod[NODE1],
1714  gNod[NODE2],
1715  gNod[NODE3], Beam::Section(iSez));
1716 
1717  /* Derivate delle derivate spaziali dei parametri di rotazione */
1718  gPrimeGrad[iSez] = InterpDeriv(gPrimeNod[NODE1],
1719  gPrimeNod[NODE2],
1720  gPrimeNod[NODE3], Beam::Section(iSez));
1721 
1722  /* Calcola le deformazioni nel sistema locale nei punti di valutazione */
1723  DefLoc[iSez] = DefLocRef[iSez]
1724  = Vec6(R[iSez].MulTV(L[iSez]) - L0[iSez],
1725  R[iSez].MulTV(Mat3x3(CGR_Rot::MatG, g[iSez])*gGrad[iSez]) + DefLocPrev[iSez].GetVec2());
1726 
1727  /* Calcola le velocita' di deformazione nel sistema locale nei punti di valutazione */
1728  DefPrimeLoc[iSez] = DefPrimeLocRef[iSez]
1729  = Vec6(R[iSez].MulTV(LPrime[iSez] + L[iSez].Cross(Omega[iSez])),
1730  R[iSez].MulTV(Mat3x3(CGR_Rot::MatG, g[iSez])*gPrimeGrad[iSez]
1731  + (Mat3x3(CGR_Rot::MatG, g[iSez])*gGrad[iSez]).Cross(Omega[iSez])));
1732 
1733  /* Calcola le azioni interne */
1734  pD[iSez]->Update(DefLoc[iSez], DefPrimeLoc[iSez]);
1735  AzLoc[iSez] = pD[iSez]->GetF();
1736 
1737  /* corregge le azioni interne locali (piezo, ecc) */
1738  AddInternalForces(AzLoc[iSez], iSez);
1739 
1740  /* Porta le azioni interne nel sistema globale */
1741  Az[iSez] = AzRef[iSez] = MultRV(AzLoc[iSez], R[iSez]);
1742 
1743  /* Aggiorna il legame costitutivo */
1744  DRef[iSez] = MultRMRt(pD[iSez]->GetFDE(), R[iSez]);
1745  ERef[iSez] = MultRMRt(pD[iSez]->GetFDEPrime(), R[iSez]);
1746  }
1747 
1748  bFirstRes = true;
1749 }
1750 
1751 doublereal
1752 ViscoElasticBeam::dGetPrivData(unsigned int i) const
1753 {
1754  ASSERT(i > 0 && i <= iGetNumPrivData());
1755 
1756  int sez = (i - 1)/iNumPrivData;
1757  int idx = (i - 1)%iNumPrivData + 1;
1758 
1759  switch (idx) {
1760  case 22:
1761  case 23:
1762  case 24:
1763  case 25:
1764  case 26:
1765  case 27:
1766  return DefPrimeLoc[sez].dGet(idx - 21);
1767 
1768  default:
1769  return Beam::dGetPrivData(i);
1770  }
1771 }
1772 
1773 /* ViscoElasticBeam - end */
1774 
1775 void
1776 ReadBeamCustomOutput(DataManager* pDM, MBDynParser& HP, unsigned int uLabel,
1777  Beam::Type BT, unsigned& uFlags, OrientationDescription& od)
1778 {
1779  uFlags = Beam::OUTPUT_NONE;
1780 
1781  while (HP.IsArg()) {
1782  unsigned uFlag;
1783 
1784  if (HP.IsKeyWord("position")) {
1785  uFlag = Beam::OUTPUT_EP_X;
1786 
1787  } else if (HP.IsKeyWord("orientation")) {
1788  uFlag = Beam::OUTPUT_EP_R;
1789 
1790  } else if (HP.IsKeyWord("configuration")) {
1791  uFlag = Beam::OUTPUT_EP_CONFIGURATION;
1792 
1793  } else if (HP.IsKeyWord("force")) {
1794  uFlag = Beam::OUTPUT_EP_F;
1795 
1796  } else if (HP.IsKeyWord("moment")) {
1797  uFlag = Beam::OUTPUT_EP_M;
1798 
1799  } else if (HP.IsKeyWord("forces")) {
1800  uFlag = Beam::OUTPUT_EP_FORCES;
1801 
1802  } else if (HP.IsKeyWord("linear" "strain")) {
1803  uFlag = Beam::OUTPUT_EP_NU;
1804 
1805  } else if (HP.IsKeyWord("angular" "strain")) {
1806  uFlag = Beam::OUTPUT_EP_K;
1807 
1808  } else if (HP.IsKeyWord("strains")) {
1809  uFlag = Beam::OUTPUT_EP_STRAINS;
1810 
1811  } else if (HP.IsKeyWord("linear" "strain" "rate")) {
1812  uFlag = Beam::OUTPUT_EP_NUP;
1813 
1814  } else if (HP.IsKeyWord("angular" "strain" "rate")) {
1815  uFlag = Beam::OUTPUT_EP_KP;
1816 
1817  } else if (HP.IsKeyWord("strain rates")) {
1818  uFlag = Beam::OUTPUT_EP_STRAIN_RATES;
1819 
1820  } else if (HP.IsKeyWord("all")) {
1821  uFlag = Beam::OUTPUT_EP_ALL;
1822 
1823  } else if (HP.IsKeyWord("none")) {
1824  uFlag = Beam::OUTPUT_NONE;
1825 
1826  } else {
1827  break;
1828  }
1829 
1830  if (uFlags & uFlag) {
1831  silent_cerr("Beam(" << uLabel << "): "
1832  "duplicate custom output "
1833  "at line " << HP.GetLineData()
1834  << std::endl);
1835  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1836  }
1837 
1838  if (uFlag & Beam::OUTPUT_EP_STRAIN_RATES) {
1839  if (BT == Beam::ELASTIC) {
1840  silent_cerr("Beam(" << uLabel << "): "
1841  "strain rates only allowed for viscoelastic beams (ignored) "
1842  "at line " << HP.GetLineData()
1843  << std::endl);
1844  uFlag &= ~Beam::OUTPUT_EP_STRAIN_RATES;
1845  }
1846  }
1847 
1848  if (uFlag & Beam::OUTPUT_EP_R) {
1849  od = ReadOptionalOrientationDescription(pDM, HP);
1850  }
1851 
1852  if (uFlag == Beam::OUTPUT_NONE) {
1853  // reset all
1854  uFlags = Beam::OUTPUT_NONE;
1855  }
1856 
1857  uFlags |= uFlag;
1858  }
1859 }
1860 
1861 
1862 void
1863 ReadOptionalBeamCustomOutput(DataManager* pDM, MBDynParser& HP, unsigned int uLabel,
1864  Beam::Type BT, unsigned& uFlags, OrientationDescription& od)
1865 {
1866  pDM->GetOutput(Elem::BEAM, uFlags, od);
1867  if (HP.IsKeyWord("custom" "output")) {
1868  ReadBeamCustomOutput(pDM, HP, uLabel, BT, uFlags, od);
1869  }
1870 }
1871 
1872 
1873 /* Legge una trave */
1874 
1875 Elem *
1876 ReadBeam(DataManager* pDM, MBDynParser& HP, unsigned int uLabel)
1877 {
1878  DEBUGCOUTFNAME("ReadBeam");
1879 
1880  /* Per ogni nodo: */
1881 
1882  /* Nodo 1 */
1883  const StructNode* pNode1 = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
1884 
1885  const Mat3x3& R1(pNode1->GetRCurr());
1886  if (HP.IsKeyWord("position")) {
1887  /* just eat it! */
1888  NO_OP;
1889  }
1890  Vec3 f1(HP.GetPosRel(ReferenceFrame(pNode1)));
1891  Mat3x3 Rn1(Eye3);
1892  if (HP.IsKeyWord("orientation") || HP.IsKeyWord("rot")) {
1893  Rn1 = HP.GetRotRel(ReferenceFrame(pNode1));
1894  }
1895 
1896  DEBUGLCOUT(MYDEBUG_INPUT, "node 1 offset (node reference frame): "
1897  << f1 << std::endl
1898  << "(global frame): "
1899  << pNode1->GetXCurr()+pNode1->GetRCurr()*f1 << std::endl);
1900 
1901  /* Nodo 2 */
1902  const StructNode* pNode2 = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
1903 
1904  Mat3x3 R2(pNode2->GetRCurr());
1905  if (HP.IsKeyWord("position")) {
1906  /* just eat it! */
1907  NO_OP;
1908  }
1909  Vec3 f2(HP.GetPosRel(ReferenceFrame(pNode2)));
1910  Mat3x3 Rn2(Eye3);
1911  if (HP.IsKeyWord("orientation") || HP.IsKeyWord("rot")) {
1912  Rn2 = HP.GetRotRel(ReferenceFrame(pNode2));
1913  }
1914 
1915  DEBUGLCOUT(MYDEBUG_INPUT, "node 2 offset (node reference frame): "
1916  << f2 << std::endl
1917  << "(global frame): "
1918  << pNode2->GetXCurr()+pNode2->GetRCurr()*f2 << std::endl);
1919 
1920  /* Nodo 3 */
1921  const StructNode* pNode3 = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
1922 
1923  Mat3x3 R3(pNode3->GetRCurr());
1924  if (HP.IsKeyWord("position")) {
1925  /* just eat it! */
1926  NO_OP;
1927  }
1928  Vec3 f3(HP.GetPosRel(ReferenceFrame(pNode3)));
1929  Mat3x3 Rn3(Eye3);
1930  if (HP.IsKeyWord("orientation") || HP.IsKeyWord("rot")) {
1931  Rn3 = HP.GetRotRel(ReferenceFrame(pNode3));
1932  }
1933 
1934  DEBUGLCOUT(MYDEBUG_INPUT, "node 3 offset (node reference frame): "
1935  << f3 << std::endl
1936  << "(global frame): "
1937  << pNode3->GetXCurr()+pNode3->GetRCurr()*f3 << std::endl);
1938 
1939  /* Per ogni punto: */
1940 
1941  /* Punto I */
1942 
1943  /* Matrice R */
1944  Mat3x3 R_I;
1945  bool b_I(false);
1946  if (HP.IsKeyWord("from" "nodes") || HP.IsKeyWord("node")) {
1947  b_I = true;
1948  } else {
1949  R_I = HP.GetRotAbs(::AbsRefFrame);
1950  }
1951 
1952  /* Legame costitutivo */
1953  ConstLawType::Type CLType_I = ConstLawType::UNKNOWN;
1954  ConstitutiveLaw6D* pD_I = HP.GetConstLaw6D(CLType_I);
1955 
1956  if (pD_I->iGetNumDof() != 0) {
1957  silent_cerr("line " << HP.GetLineData()
1958  << ": Beam(" << uLabel << ") "
1959  "does not support dynamic constitutive laws yet"
1960  << std::endl);
1961  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1962  }
1963 
1964 #ifdef DEBUG
1965  Mat6x6 MTmp(pD_I->GetFDE());
1966  Mat3x3 D11(MTmp.GetMat11());
1967  Mat3x3 D12(MTmp.GetMat12());
1968  Mat3x3 D21(MTmp.GetMat21());
1969  Mat3x3 D22(MTmp.GetMat22());
1970 
1971  DEBUGLCOUT(MYDEBUG_INPUT,
1972  "First point matrix D11: " << std::endl << D11 << std::endl
1973  << "First point matrix D12: " << std::endl << D12 << std::endl
1974  << "First point matrix D21: " << std::endl << D21 << std::endl
1975  << "First point matrix D22: " << std::endl << D22 << std::endl);
1976 #endif
1977 
1978 
1979  /* Punto II */
1980 
1981  /* Matrice R */
1982  Mat3x3 RII;
1983  bool bII(false);
1984  if (HP.IsKeyWord("same")) {
1985  if (b_I) {
1986  bII = true;
1987  } else {
1988  RII = R_I;
1989  }
1990  } else {
1991  if (HP.IsKeyWord("from" "nodes") || HP.IsKeyWord("node")) {
1992  bII = true;
1993  } else {
1994  RII = HP.GetRotAbs(::AbsRefFrame);
1995  }
1996  }
1997 
1998  /* Legame costitutivo */
1999  ConstLawType::Type CLTypeII = ConstLawType::UNKNOWN;
2000  ConstitutiveLaw6D* pDII = NULL;
2001 
2002  if (HP.IsKeyWord("same")) {
2003  pDII = pD_I->pCopy();
2004  CLTypeII = CLType_I;
2005 
2006  } else {
2007  pDII = HP.GetConstLaw6D(CLTypeII);
2008 
2009  if (pDII->iGetNumDof() != 0) {
2010  silent_cerr("line " << HP.GetLineData()
2011  << ": Beam(" << uLabel << ") "
2012  "does not support dynamic constitutive laws yet"
2013  << std::endl);
2014  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2015  }
2016  }
2017 
2018  Beam::Type Type;
2019  if (CLType_I == ConstLawType::ELASTIC && CLTypeII == ConstLawType::ELASTIC) {
2020  Type = Beam::ELASTIC;
2021  } else {
2022  Type = Beam::VISCOELASTIC;
2023  }
2024 
2025 #ifdef DEBUG
2026  MTmp = pDII->GetFDE();
2027  D11 = MTmp.GetMat11();
2028  D12 = MTmp.GetMat12();
2029  D21 = MTmp.GetMat21();
2030  D22 = MTmp.GetMat22();
2031 
2032  DEBUGLCOUT(MYDEBUG_INPUT,
2033  "Second point matrix D11: " << std::endl << D11 << std::endl
2034  << "Second point matrix D12: " << std::endl << D12 << std::endl
2035  << "Second point matrix D21: " << std::endl << D21 << std::endl
2036  << "Second point matrix D22: " << std::endl << D22 << std::endl);
2037 #endif
2038 
2039  flag fPiezo(0);
2040  Mat3xN PiezoMat[2][2];
2041  integer iNumElec = 0;
2042  const ScalarDifferentialNode **pvElecDofs = 0;
2043  if (HP.IsKeyWord("piezoelectric" "actuator")) {
2044  fPiezo = flag(1);
2045  DEBUGLCOUT(MYDEBUG_INPUT,
2046  "Piezoelectric actuator beam is expected" << std::endl);
2047 
2048  iNumElec = HP.GetInt();
2049  DEBUGLCOUT(MYDEBUG_INPUT,
2050  "piezo actuator " << uLabel << " has " << iNumElec
2051  << " electrodes" << std::endl);
2052  if (iNumElec <= 0) {
2053  silent_cerr("PiezoElectricBeam(" << uLabel << "): "
2054  "illegal number of electric nodes " << iNumElec
2055  << " at line " << HP.GetLineData() << std::endl);
2056  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2057  }
2058 
2059  SAFENEWARR(pvElecDofs, const ScalarDifferentialNode *, iNumElec);
2060 
2061  for (integer i = 0; i < iNumElec; i++) {
2062  pvElecDofs[i] = pDM->ReadNode<const ScalarDifferentialNode, Node::ABSTRACT>(HP);
2063  }
2064 
2065  PiezoMat[0][0].Resize(iNumElec);
2066  PiezoMat[1][0].Resize(iNumElec);
2067  PiezoMat[0][1].Resize(iNumElec);
2068  PiezoMat[1][1].Resize(iNumElec);
2069 
2070  /* leggere le matrici (6xN sez. 1, 6xN sez. 2) */
2071  HP.GetMat6xN(PiezoMat[0][0], PiezoMat[1][0], iNumElec);
2072  if (HP.IsKeyWord("same")) {
2073  PiezoMat[0][1].Copy(PiezoMat[0][0]);
2074  PiezoMat[1][1].Copy(PiezoMat[1][0]);
2075  } else {
2076  HP.GetMat6xN(PiezoMat[0][1], PiezoMat[1][1], iNumElec);
2077  }
2078 
2079 #if 0
2080  DEBUGLCOUT(MYDEBUG_INPUT, "Piezo matrix I:" << std::endl << PiezoMat[0][0] << PiezoMat[1][0]);
2081  DEBUGLCOUT(MYDEBUG_INPUT, "Piezo matrix II:" << std::endl << PiezoMat[0][1] << PiezoMat[1][1]);
2082 #endif /* 0 */
2083  }
2084 
2085  OrientationDescription od = UNKNOWN_ORIENTATION_DESCRIPTION;
2086  unsigned uFlags = Beam::OUTPUT_NONE;
2087  ReadOptionalBeamCustomOutput(pDM, HP, uLabel, Type, uFlags, od);
2088 
2089  flag fOut = pDM->fReadOutput(HP, Elem::BEAM);
2090  if (fOut) {
2091  fOut |= uFlags;
2092  }
2093 
2094  /* Se necessario, interpola i parametri di rotazione delle sezioni */
2095  if (b_I || bII) {
2096  Mat3x3 R(R2*Rn2);
2097  Vec3 g1(CGR_Rot::Param, R.MulTM(R1*Rn1));
2098  Vec3 g3(CGR_Rot::Param, R.MulTM(R3*Rn3));
2099  if (b_I) {
2100  R_I = R*Mat3x3(CGR_Rot::MatR, Beam::InterpState(g1, Zero3, g3, Beam::S_I));
2101  }
2102  if (bII) {
2103  RII = R*Mat3x3(CGR_Rot::MatR, Beam::InterpState(g1, Zero3, g3, Beam::SII));
2104  }
2105  }
2106 
2107  std::ostream& out = pDM->GetLogFile();
2108  out << "beam3: " << uLabel
2109  << " " << pNode1->GetLabel()
2110  << " ", f1.Write(out, " ")
2111  << " " << pNode2->GetLabel()
2112  << " ", f2.Write(out, " ")
2113  << " " << pNode3->GetLabel()
2114  << " ", f3.Write(out, " ")
2115  << std::endl;
2116 
2117  Elem* pEl = NULL;
2118 
2119  if ((CLType_I == ConstLawType::ELASTIC)
2120  && (CLTypeII == ConstLawType::ELASTIC))
2121  {
2122  if (fPiezo == 0) {
2123  SAFENEWWITHCONSTRUCTOR(pEl,
2124  Beam,
2125  Beam(uLabel,
2126  pNode1, pNode2, pNode3,
2127  f1, f2, f3,
2128  Rn1, Rn2, Rn3,
2129  R_I, RII,
2130  pD_I, pDII,
2131  od, fOut));
2132  } else {
2133  SAFENEWWITHCONSTRUCTOR(pEl,
2134  PiezoActuatorBeam,
2135  PiezoActuatorBeam(uLabel,
2136  pNode1, pNode2, pNode3,
2137  f1, f2, f3,
2138  Rn1, Rn2, Rn3,
2139  R_I, RII,
2140  pD_I, pDII,
2141  iNumElec,
2142  pvElecDofs,
2143  PiezoMat[0][0], PiezoMat[1][0],
2144  PiezoMat[0][1], PiezoMat[1][1],
2145  od, fOut));
2146  }
2147 
2148  } else /* At least one is VISCOUS or VISCOELASTIC */ {
2149  if (fPiezo == 0) {
2150  SAFENEWWITHCONSTRUCTOR(pEl,
2151  ViscoElasticBeam,
2152  ViscoElasticBeam(uLabel,
2153  pNode1, pNode2, pNode3,
2154  f1, f2, f3,
2155  Rn1, Rn2, Rn3,
2156  R_I, RII,
2157  pD_I, pDII,
2158  od, fOut));
2159  } else {
2160  SAFENEWWITHCONSTRUCTOR(pEl,
2161  PiezoActuatorVEBeam,
2162  PiezoActuatorVEBeam(uLabel,
2163  pNode1, pNode2, pNode3,
2164  f1, f2, f3,
2165  Rn1, Rn2, Rn3,
2166  R_I, RII,
2167  pD_I, pDII,
2168  iNumElec,
2169  pvElecDofs,
2170  PiezoMat[0][0], PiezoMat[1][0],
2171  PiezoMat[0][1], PiezoMat[1][1],
2172  od, fOut));
2173  }
2174  }
2175 
2176  /* Costruttore normale
2177  * Beam(unsigned int uL,
2178  * const StructNode* pN1, const StructNode* pN2, const StructNode* pN3,
2179  * const Vec3& X1, const Vec3& X2, const Vec3& X3,
2180  * const Vec3& F1, const Vec3& F2, const Vec3& F3,
2181  * const Mat3x3& r_I, const Mat3x3& rII,
2182  * const Mat3x3& d11_I, const Mat3x3& d12_I,
2183  * const Mat3x3& d21_I, const Mat3x3& d22_I,
2184  * const Mat3x3& d11II, const Mat3x3& d12II,
2185  * const Mat3x3& d21II, const Mat3x3& d22II,
2186  * const Vec3& eps0_I, const Vec3& k0_I,
2187  * const Vec3& eps0II, const Vec3& k0II);
2188  */
2189 
2190 
2191  /* Se non c'e' il punto e virgola finale */
2192  if (HP.IsArg()) {
2193  silent_cerr("semicolon expected "
2194  "at line " << HP.GetLineData() << std::endl);
2195  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2196  }
2197 
2198  return pEl;
2199 } /* End of ReadBeam() */
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
CGR_Rot::MatGm1
const MatGm1_Manip MatGm1
Definition: matvec3.cc:647
Beam::Beam
Beam(unsigned int uL, const StructNode *pN1, const StructNode *pN2, const StructNode *pN3, const Vec3 &F1, const Vec3 &F2, const Vec3 &F3, const Mat3x3 &R1, const Mat3x3 &R2, const Mat3x3 &R3, const Mat3x3 &r_I, const Mat3x3 &rII, const ConstitutiveLaw6D *pD_I, const ConstitutiveLaw6D *pDII, OrientationDescription ood, flag fOut)
Definition: beam.cc:64
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
ViscoElasticBeam::AssStiffnessVec
virtual void AssStiffnessVec(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:1522
MultRMRt
Mat3x3 MultRMRt(const Mat3x3 &m, const Mat3x3 &R)
Definition: matvec3.cc:1162
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
Beam::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:886
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
ViscoElasticBeam::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: beam.cc:1384
Vec3::Write
std::ostream & Write(std::ostream &out, const char *sFill=" ") const
Definition: matvec3.cc:738
MultRV
Vec3 MultRV(const Vec3 &v, const Mat3x3 &R)
Definition: matvec3.cc:1144
Beam::iGetNumPrivData
virtual unsigned int iGetNumPrivData(void) const
Definition: beam.cc:263
flag
long int flag
Definition: mbdyn.h:43
StructNode::GetRRef
virtual const Mat3x3 & GetRRef(void) const
Definition: strnode.h:1006
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
Beam::AssInertiaMat
virtual void AssInertiaMat(FullSubMatrixHandler &, FullSubMatrixHandler &, doublereal, const VectorHandler &, const VectorHandler &)
Definition: beam.h:232
Beam::DefLocPrev
Vec6 DefLocPrev[NUMSEZ]
Definition: beam.h:182
ConstitutiveLawOwner::pGetConstLaw
ConstitutiveLaw< T, Tder > * pGetConstLaw(void) const
Definition: constltp.h:278
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
Beam::Init
void Init(void)
Definition: beam.cc:191
DEBUGCOUTFNAME
#define DEBUGCOUTFNAME(fname)
Definition: myassert.h:256
VectorHandler::ResizeReset
virtual void ResizeReset(integer)
Definition: vh.cc:55
HighParser::GetInt
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
Vec6::GetVec2
const Vec3 & GetVec2(void) const
Definition: matvec6.h:76
Beam::f
const Vec3 f[NUMNODES]
Definition: beam.h:145
MatCross
const MatCross_Manip MatCross
Definition: matvec3.cc:639
OutputHandler::UseNetCDF
bool UseNetCDF(int out) const
Definition: output.cc:491
Beam::InterpState
static Vec3 InterpState(const Vec3 &v1, const Vec3 &v2, const Vec3 &v3, enum Section Sec)
Definition: beam.cc:429
StructNode::GetgCurr
virtual const Vec3 & GetgCurr(void) const
Definition: strnode.h:982
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
MBDynParser::GetRotAbs
Mat3x3 GetRotAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1857
StructDispNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: strnode.cc:145
Beam::InitialAssRes
virtual SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: beam.cc:1265
Beam::L
Vec3 L[NUMSEZ]
Definition: beam.h:190
Vec3::Dot
doublereal Dot(const Vec3 &v) const
Definition: matvec3.h:243
FullSubMatrixHandler::Add
void Add(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:209
Beam::AddInternalForces
virtual void AddInternalForces(Vec6 &, unsigned int)
Definition: beam.h:227
ConstitutiveLaw6DOwner
ConstitutiveLawOwner< Vec6, Mat6x6 > ConstitutiveLaw6DOwner
Definition: constltp.h:380
Write
std::ostream & Write(std::ostream &out, const FullMatrixHandler &m, const char *s, const char *s2)
Definition: fullmh.cc:376
Beam::RNode
const Mat3x3 RNode[NUMNODES]
Definition: beam.h:147
Beam::SetValue
void SetValue(DataManager *pDM, VectorHandler &, VectorHandler &, SimulationEntity::Hints *ph=0)
Definition: beam.cc:966
Beam::AssInertiaVec
virtual void AssInertiaVec(SubVectorHandler &, doublereal, const VectorHandler &, const VectorHandler &)
Definition: beam.h:241
MBDynParser::GetMat6xN
virtual void GetMat6xN(Mat3xN &m1, Mat3xN &m2, integer iNumCols)
Definition: mbpar.cc:2743
ViscoElasticBeam::ViscoElasticBeam
ViscoElasticBeam(unsigned int uL, const StructNode *pN1, const StructNode *pN2, const StructNode *pN3, const Vec3 &F1, const Vec3 &F2, const Vec3 &F3, const Mat3x3 &R1, const Mat3x3 &R2, const Mat3x3 &R3, const Mat3x3 &r_I, const Mat3x3 &rII, const ConstitutiveLaw6D *pD_I, const ConstitutiveLaw6D *pDII, OrientationDescription ood, flag fOut)
Definition: beam.cc:1310
Beam::Omega
Vec3 Omega[NUMSEZ]
Definition: beam.h:173
Beam::pD
ConstitutiveLaw6DOwner * pD[NUMSEZ]
Definition: beam.h:156
Beam::AssStiffnessVec
virtual void AssStiffnessVec(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:697
ConstitutiveLaw::GetFDE
virtual const Tder & GetFDE(void) const
Definition: constltp.h:129
VectorHandler::Sub
virtual void Sub(integer iRow, const Vec3 &v)
Definition: vh.cc:78
Eye3
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
OrientationDescription
OrientationDescription
Definition: matvec3.h:1597
iNumPrivData
static const unsigned int iNumPrivData
Definition: beam.cc:249
DataManager::GetOutput
void GetOutput(Elem::Type t, unsigned &, OrientationDescription &) const
Definition: dataman.cc:871
StructNode::GetWRef
virtual const Vec3 & GetWRef(void) const
Definition: strnode.h:1024
mb_deye< Mat3x3 >
Mat3x3 mb_deye< Mat3x3 >(const doublereal d)
Definition: matvec3.h:1584
NO_OP
#define NO_OP
Definition: myassert.h:74
ReadOptionalOrientationDescription
OrientationDescription ReadOptionalOrientationDescription(DataManager *pDM, MBDynParser &HP)
Definition: dataman3.cc:2531
Beam::DefPrimeLoc
Vec6 DefPrimeLoc[NUMSEZ]
Definition: beam.h:185
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
Beam::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: beam.cc:595
Beam::AfterPredict
virtual void AfterPredict(VectorHandler &, VectorHandler &)
Definition: beam.cc:988
Beam::iNumPrivData
static const unsigned int iNumPrivData
Definition: beam.h:100
ViscoElasticBeam::DefPrimeLocRef
Vec6 DefPrimeLocRef[NUMSEZ]
Definition: beam.h:461
Mat3x3::GetVec
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
Beam::pGetNode
virtual const StructNode * pGetNode(unsigned int i) const
Definition: beam.cc:1290
Beam::RPrev
Mat3x3 RPrev[NUMSEZ]
Definition: beam.h:153
Beam::iGetPrivDataIdx
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: beam.cc:346
RotManip::VecRot
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
ViscoElasticBeam::dGetPrivData
virtual doublereal dGetPrivData(unsigned int i) const
Definition: beam.cc:1752
Beam::DefLocRef
Vec6 DefLocRef[NUMSEZ]
Definition: beam.h:181
ReadBeamCustomOutput
void ReadBeamCustomOutput(DataManager *pDM, MBDynParser &HP, unsigned int uLabel, Beam::Type BT, unsigned &uFlags, OrientationDescription &od)
Definition: beam.cc:1776
Mat3x3::MulTV
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
Mat6x6::GetMat11
Mat3x3 GetMat11(void)
Definition: matvec6.h:320
Beam::Output
virtual void Output(OutputHandler &OH) const
Definition: beam.cc:1141
Beam::od
OrientationDescription od
Definition: beam.h:116
ViscoElasticBeam::AssStiffnessMat
virtual void AssStiffnessMat(FullSubMatrixHandler &WMA, FullSubMatrixHandler &WMB, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:1404
SubVectorHandler::PutRowIndex
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
Zero3x3
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
Beam::AzLoc
Vec6 AzLoc[NUMSEZ]
Definition: beam.h:179
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)
ViscoElasticBeam::SetValue
void SetValue(DataManager *pDM, VectorHandler &, VectorHandler &, SimulationEntity::Hints *ph=0)
Definition: beam.cc:1635
Vec6
Definition: matvec6.h:37
Beam::LRef
Vec3 LRef[NUMSEZ]
Definition: beam.h:192
Beam::L0
Vec3 L0[NUMSEZ]
Definition: beam.h:189
Vec6::GetVec1
const Vec3 & GetVec1(void) const
Definition: matvec6.h:72
ViscoElasticBeam::ERef
Mat6x6 ERef[NUMSEZ]
Definition: beam.h:463
MatR2EulerAngles313
Vec3 MatR2EulerAngles313(const Mat3x3 &R)
Definition: matvec3.cc:927
Mat6x6::SubMat12
void SubMat12(const Mat3x3 &x)
Definition: matvec6.h:416
MBDynParser::GetPosRel
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
Beam::OmegaRef
Vec3 OmegaRef[NUMSEZ]
Definition: beam.h:174
Mat6x6::SubMat22
void SubMat22(const Mat3x3 &x)
Definition: matvec6.h:420
Beam::~Beam
virtual ~Beam(void)
Definition: beam.cc:236
ViscoElasticBeam::Init
void Init(void)
Definition: beam.cc:1363
VariableSubMatrixHandler::SetNullMatrix
void SetNullMatrix(void)
Definition: submat.h:1159
dN3
const doublereal dN3[2][3]
Definition: shapefnc.cc:157
false
Definition: mbdyn.h:76
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
ConstitutiveLaw::pCopy
virtual ConstitutiveLaw< T, Tder > * pCopy(void) const =0
Zero6
const Vec6 Zero6(0., 0., 0., 0., 0., 0.)
StructDispNode::od
OrientationDescription od
Definition: strnode.h:111
CGR_Rot::MatG
const MatG_Manip MatG
Definition: matvec3.cc:646
OutputHandler::GetCurrentStep
long GetCurrentStep(void) const
Definition: output.h:116
Vec3::dGet
const doublereal & dGet(unsigned short int iRow) const
Definition: matvec3.h:285
Beam::R
Mat3x3 R[NUMSEZ]
Definition: beam.h:151
Beam
Definition: beam.h:56
MatR2EulerAngles123
Vec3 MatR2EulerAngles123(const Mat3x3 &R)
Definition: matvec3.cc:893
true
Definition: mbdyn.h:77
Beam::dGetPrivData
virtual doublereal dGetPrivData(unsigned int i) const
Definition: beam.cc:381
ViscoElasticBeam::LPrime
Vec3 LPrime[NUMSEZ]
Definition: beam.h:453
StructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const =0
StructDispNode::iGetFirstPositionIndex
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
StructNode::GetWCurr
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
Beam::Section
Section
Definition: beam.h:136
dRaDegr
const doublereal dRaDegr
Definition: matvec3.cc:884
WithLabel::uLabel
unsigned int uLabel
Definition: withlab.h:44
ViscoElasticBeam::AfterPredict
virtual void AfterPredict(VectorHandler &, VectorHandler &)
Definition: beam.cc:1657
Beam::InterpDeriv
Vec3 InterpDeriv(const Vec3 &v1, const Vec3 &v2, const Vec3 &v3, enum Section Sec)
Definition: beam.cc:446
OutputHandler::IsOpen
bool IsOpen(int out) const
Definition: output.cc:395
Beam::Type
Type
Definition: beam.h:64
ConstitutiveLawOwner::AfterConvergence
void AfterConvergence(const T &Eps, const T &EpsPrime=mb_zero< T >())
Definition: constltp.h:288
Beam::Omega0
virtual void Omega0(void)
Definition: beam.cc:519
CGR_Rot::MatR
const MatR_Manip MatR
Definition: matvec3.cc:645
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
Beam::pNode
const StructNode * pNode[NUMNODES]
Definition: beam.h:142
MBDynParser::GetConstLaw6D
ConstitutiveLaw6D * GetConstLaw6D(ConstLawType::Type &clt)
Definition: mbpar.cc:1995
Beam::iGetPrivDataIdx_int
static unsigned int iGetPrivDataIdx_int(const char *s, ConstLawType::Type type)
Definition: beam.cc:269
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
Mat3x3::MulTM
Mat3x3 MulTM(const Mat3x3 &m) const
Definition: matvec3.cc:500
grad::Cross
VectorExpression< VectorCrossExpr< VectorLhsExpr, VectorRhsExpr >, 3 > Cross(const VectorExpression< VectorLhsExpr, 3 > &u, const VectorExpression< VectorRhsExpr, 3 > &v)
Definition: matvec.h:3248
ViscoElasticBeam::gPrime
Vec3 gPrime[NUMSEZ]
Definition: beam.h:454
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
Beam::DsDxi
virtual void DsDxi(void)
Definition: beam.cc:466
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
VectorHandler::Add
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
Mat3xN::Resize
void Resize(integer ns)
Definition: matvec3n.cc:318
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
ReadOptionalBeamCustomOutput
void ReadOptionalBeamCustomOutput(DataManager *pDM, MBDynParser &HP, unsigned int uLabel, Beam::Type BT, unsigned &uFlags, OrientationDescription &od)
Definition: beam.cc:1863
ConstitutiveLaw::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: constltp.h:138
Beam::bFirstRes
bool bFirstRes
Definition: beam.h:197
Beam::AzRef
Vec6 AzRef[NUMSEZ]
Definition: beam.h:178
MatCrossCross
const MatCrossCross_Manip MatCrossCross
Definition: matvec3.cc:640
Beam::RRef
Mat3x3 RRef[NUMSEZ]
Definition: beam.h:152
Beam::g
Vec3 g[NUMSEZ]
Definition: beam.h:188
OutputHandler::Beams
std::ostream & Beams(void) const
Definition: output.h:457
STRLENOF
#define STRLENOF(s)
Definition: mbdyn.h:166
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
Elem
Definition: elem.h:75
Beam::AssMats
void AssMats(VariableSubMatrixHandler &WorkMatA, VariableSubMatrixHandler &WorkMatB, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:822
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
ConstitutiveLawOwner::GetF
const T & GetF(void) const
Definition: constltp.h:293
Vec3::pGetVec
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
StructNode::GetgRef
virtual const Vec3 & GetgRef(void) const
Definition: strnode.h:976
Beam::InitialAssJac
virtual VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
Definition: beam.cc:1236
DataManager::GetLogFile
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
Beam::bInverseDynamics
virtual bool bInverseDynamics(void) const
Definition: beam.cc:959
ToBeOutput::fToBeOutput
virtual flag fToBeOutput(void) const
Definition: output.cc:884
MatR2EulerAngles321
Vec3 MatR2EulerAngles321(const Mat3x3 &R)
Definition: matvec3.cc:941
StructNode::GetgPCurr
virtual const Vec3 & GetgPCurr(void) const
Definition: strnode.h:994
StructDispNode::XCurr
Vec3 XCurr
Definition: strnode.h:90
StructDispNode::GetVCurr
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
FullSubMatrixHandler::Sub
void Sub(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:215
SAFENEWARR
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
Beam::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: beam.cc:570
ReadBeam
Elem * ReadBeam(DataManager *pDM, MBDynParser &HP, unsigned int uLabel)
Definition: beam.cc:1876
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
Beam::OutputPrepare
virtual void OutputPrepare(OutputHandler &OH)
Definition: beam.cc:1047
StructNode::GetgPRef
virtual const Vec3 & GetgPRef(void) const
Definition: strnode.h:988
Beam::Az
Vec6 Az[NUMSEZ]
Definition: beam.h:177
Mat3xN::Copy
const Mat3xN & Copy(const Mat3xN &m)
Definition: matvec3n.cc:407
Beam::AssJac
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:772
doublereal
double doublereal
Definition: colamd.c:52
Beam::bConsistentInertia
const bool bConsistentInertia
Definition: beam.h:162
integer
long int integer
Definition: colamd.c:51
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
ConstitutiveLawOwner::Update
void Update(const T &Eps, const T &EpsPrime=mb_zero< T >())
Definition: constltp.h:283
Beam::AssStiffnessMat
virtual void AssStiffnessMat(FullSubMatrixHandler &WMA, FullSubMatrixHandler &WMB, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beam.cc:613
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
DataManager::ReadNode
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
Vec6::dGet
const doublereal & dGet(unsigned short int i) const
Definition: matvec6.h:182
DEBUGLCOUT
#define DEBUGLCOUT(level, msg)
Definition: myassert.h:244
ViscoElasticBeam::LPrimeRef
Vec3 LPrimeRef[NUMSEZ]
Definition: beam.h:456
Beam::DRef
Mat6x6 DRef[NUMSEZ]
Definition: beam.h:159
dN3P
const doublereal dN3P[2][3]
Definition: shapefnc.cc:162
Beam::p
Vec3 p[NUMSEZ]
Definition: beam.h:187
R
Mat3x3 R
Definition: mbdyn_util_oct.cc:74
OutputHandler::UseText
bool UseText(int out) const
Definition: output.cc:446
Beam::DefLoc
Vec6 DefLoc[NUMSEZ]
Definition: beam.h:180
Beam::Restart_
virtual std::ostream & Restart_(std::ostream &out) const
Definition: beam.cc:576
Beam::GetBeamType
virtual Beam::Type GetBeamType(void) const
Definition: beam.h:292
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
Beam::dsdxi
doublereal dsdxi[NUMSEZ]
Definition: beam.h:194