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 */
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)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 */
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,doublereal dM_I,const Vec3 & s0_I,const Mat3x3 & j0_I,doublereal dMII,const Vec3 & s0II,const Mat3x3 & j0II,OrientationDescription ood,flag fOut)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
Init(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
~Beam(void)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
iGetNumPrivData(void) const263 Beam::iGetNumPrivData(void) const
264 {
265 return 2*iNumPrivData;
266 }
267
268 unsigned int
iGetPrivDataIdx_int(const char * s,ConstLawType::Type type)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
iGetPrivDataIdx(const char * s) const346 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
dGetPrivData(unsigned int i) const381 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
InterpState(const Vec3 & v1,const Vec3 & v2,const Vec3 & v3,enum Section Sec)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
InterpDeriv(const Vec3 & v1,const Vec3 & v2,const Vec3 & v3,enum Section Sec)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
DsDxi(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
Omega0(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() */
2200