1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aeroelem.cc,v 1.124 2017/01/12 14:45:58 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31
32 /* Elementi aerodinamici bidimensionali */
33
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35
36 #include "dataman.h"
37 #include "aerodata.h"
38 #include "aeroelem.h"
39
40 #include "shapefnc.h"
41
42 #include "aerodc81.h"
43 #include "c81data.h"
44 #include "Rot.hh"
45
46 #include <sstream>
47
48 /* AerodynamicOutput - begin */
49
AerodynamicOutput(flag f,int iNP,OrientationDescription ood)50 AerodynamicOutput::AerodynamicOutput(flag f, int iNP,
51 OrientationDescription ood)
52 : m_eOutput(f),
53 od(ood)
54 {
55 SetOutputFlag(f, iNP);
56 }
57
~AerodynamicOutput(void)58 AerodynamicOutput::~AerodynamicOutput(void)
59 {
60 NO_OP;
61 }
62
63 void
SetOutputFlag(flag f,int iNP)64 AerodynamicOutput::SetOutputFlag(flag f, int iNP)
65 {
66 m_eOutput = f;
67
68 if (IsOutput() && IsPGAUSS()) {
69 OutputData.resize(iNP);
70
71 } else {
72 OutputData.resize(0);
73 }
74 }
75
76 void
ResetIterator(void)77 AerodynamicOutput::ResetIterator(void)
78 {
79 if (IsOutput() && IsPGAUSS()) {
80 ASSERT(!OutputData.empty());
81 OutputIter = OutputData.begin();
82 }
83
84 #ifdef USE_NETCDF
85 if (!NetCDFOutputData.empty()) {
86 NetCDFOutputIter = NetCDFOutputData.begin();
87 }
88 #endif // USE_NETCDF
89 }
90
91 void
SetData(const Vec3 & v,const doublereal * pd,const Vec3 & X,const Mat3x3 & R,const Vec3 & V,const Vec3 & W,const Vec3 & F,const Vec3 & M)92 AerodynamicOutput::SetData(const Vec3& v, const doublereal* pd,
93 const Vec3& X, const Mat3x3& R, const Vec3& V, const Vec3& W, const Vec3& F, const Vec3& M)
94 {
95 if (IsPGAUSS()) {
96 ASSERT(!OutputData.empty());
97 ASSERT(OutputIter >= OutputData.begin());
98 ASSERT(OutputIter < OutputData.end());
99
100 OutputIter->alpha = 180./M_PI*atan2(-v(2), v(1));
101 OutputIter->f = Vec3(pd[1], pd[0], pd[5]);
102
103 // move iterator forward
104 ++OutputIter;
105 }
106
107 #ifdef USE_NETCDF
108 if (!NetCDFOutputData.empty()) {
109 ASSERT(NetCDFOutputIter >= NetCDFOutputData.begin());
110 ASSERT(NetCDFOutputIter < NetCDFOutputData.end());
111
112 if (NetCDFOutputIter->Var_X) NetCDFOutputIter->X = X;
113 if (NetCDFOutputIter->Var_Phi) NetCDFOutputIter->R = R;
114 if (NetCDFOutputIter->Var_V) NetCDFOutputIter->V = V;
115 if (NetCDFOutputIter->Var_W) NetCDFOutputIter->W = W;
116 if (NetCDFOutputIter->Var_F) NetCDFOutputIter->F = F;
117 if (NetCDFOutputIter->Var_M) NetCDFOutputIter->M = M;
118
119 ++NetCDFOutputIter;
120 }
121 #endif // USE_NETCDF
122 }
123
124 AerodynamicOutput::eOutput
GetOutput(void) const125 AerodynamicOutput::GetOutput(void) const
126 {
127 return eOutput(m_eOutput & AEROD_OUT_MASK);
128 }
129
130 bool
IsOutput(void) const131 AerodynamicOutput::IsOutput(void) const
132 {
133 return (m_eOutput & ToBeOutput::OUTPUT);
134 }
135
136 bool
IsSTD(void) const137 AerodynamicOutput::IsSTD(void) const
138 {
139 return GetOutput() == AEROD_OUT_STD;
140 }
141
142 bool
IsPGAUSS(void) const143 AerodynamicOutput::IsPGAUSS(void) const
144 {
145 return GetOutput() == AEROD_OUT_PGAUSS;
146 }
147
148 bool
IsNODE(void) const149 AerodynamicOutput::IsNODE(void) const
150 {
151 return GetOutput() == AEROD_OUT_NODE;
152 }
153
154 /* AerodynamicOutput - end */
155
156
157 /* Aerodynamic2DElem - begin */
158
159 static const bool bDefaultUseJacobian = false;
160
161 template <unsigned iNN>
Aerodynamic2DElem(unsigned int uLabel,const DofOwner * pDO,InducedVelocity * pR,bool bPassive,const Shape * pC,const Shape * pF,const Shape * pV,const Shape * pT,const Shape * pTL,integer iN,AeroData * a,const DriveCaller * pDC,bool bUseJacobian,OrientationDescription ood,flag fOut)162 Aerodynamic2DElem<iNN>::Aerodynamic2DElem(unsigned int uLabel,
163 const DofOwner *pDO,
164 InducedVelocity* pR, bool bPassive,
165 const Shape* pC, const Shape* pF,
166 const Shape* pV, const Shape* pT,
167 const Shape* pTL,
168 integer iN, AeroData* a,
169 const DriveCaller* pDC,
170 bool bUseJacobian,
171 OrientationDescription ood,
172 flag fOut)
173 : Elem(uLabel, fOut),
174 AerodynamicElem(uLabel, pDO, fOut),
175 InitialAssemblyElem(uLabel, fOut),
176 DriveOwner(pDC),
177 AerodynamicOutput(fOut, iNN*iN, ood),
178 aerodata(a),
179 pIndVel(pR),
180 bPassiveInducedVelocity(bPassive),
181 Chord(pC),
182 ForcePoint(pF),
183 VelocityPoint(pV),
184 Twist(pT),
185 TipLoss(pTL),
186 GDI(iN),
187 OUTA(iNN*iN, outa_Zero),
188 bJacobian(bUseJacobian)
189 {
190 DEBUGCOUTFNAME("Aerodynamic2DElem::Aerodynamic2DElem");
191
192 ASSERT(iNN >= 1 && iNN <= 3);
193 ASSERT(aerodata != 0);
194 }
195
196 template <unsigned iNN>
~Aerodynamic2DElem(void)197 Aerodynamic2DElem<iNN>::~Aerodynamic2DElem(void)
198 {
199 DEBUGCOUTFNAME("Aerodynamic2DElem::~Aerodynamic2DElem");
200
201 SAFEDELETE(aerodata);
202 }
203
204 /*
205 * overload della funzione di ToBeOutput();
206 * serve per allocare il vettore dei dati di output se il flag
207 * viene settato dopo la costruzione
208 */
209 template <unsigned iNN>
210 void
SetOutputFlag(flag f)211 Aerodynamic2DElem<iNN>::SetOutputFlag(flag f)
212 {
213 DEBUGCOUTFNAME("Aerodynamic2DElem::SetOutputFlag");
214 ToBeOutput::SetOutputFlag(f);
215 AerodynamicOutput::SetOutputFlag(f, iNN*GDI.iGetNum());
216 }
217
218 /* Dimensioni del workspace */
219 template <unsigned iNN>
220 void
WorkSpaceDim(integer * piNumRows,integer * piNumCols) const221 Aerodynamic2DElem<iNN>::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
222 {
223 *piNumRows = *piNumCols = iNN*6 + iGetNumDof();
224 }
225
226 /* inherited from SimulationEntity */
227 template <unsigned iNN>
228 unsigned int
iGetNumDof(void) const229 Aerodynamic2DElem<iNN>::iGetNumDof(void) const
230 {
231 return aerodata->iGetNumDof()*iNN*GDI.iGetNum();
232 }
233
234 template <unsigned iNN>
235 std::ostream&
DescribeDof(std::ostream & out,const char * prefix,bool bInitial) const236 Aerodynamic2DElem<iNN>::DescribeDof(std::ostream& out,
237 const char *prefix, bool bInitial) const
238 {
239 ASSERT(bInitial == false);
240
241 integer iFirstIndex = iGetFirstIndex();
242 integer iNumDof = aerodata->iGetNumDof();
243
244 for (unsigned b = 0; b < iNN; b++) {
245 for (integer i = 0; i < GDI.iGetNum(); i++) {
246 out
247 << prefix << iFirstIndex + 1 << "->" << iFirstIndex + iNumDof << ": "
248 "internal states at point #" << i << " of " << GDI.iGetNum() << ", "
249 "block #" << b << " of " << iNN << std::endl;
250
251 iFirstIndex += iNumDof;
252 }
253 }
254
255 return out;
256 }
257
258 static const char *elemnames[] = {
259 "AerodynamicBody",
260 "AerodynamicBeam2",
261 "AerodynamicBeam3",
262 0
263 };
264
265 template <unsigned iNN>
266 void
DescribeDof(std::vector<std::string> & desc,bool bInitial,int i) const267 Aerodynamic2DElem<iNN>::DescribeDof(std::vector<std::string>& desc,
268 bool bInitial, int i) const
269 {
270 if (i < -1) {
271 // error
272 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
273 }
274
275 ASSERT(bInitial == false);
276
277 std::ostringstream os;
278 os << elemnames[iNN - 1] << "(" << GetLabel() << ")";
279
280 integer iNumDof = aerodata->iGetNumDof();
281 if (i == -1) {
282 integer iTotDof = iNumDof*GDI.iGetNum();
283 desc.resize(iTotDof);
284
285 std::string name(os.str());
286
287 for (integer s = 0; s < iTotDof; s++) {
288 os.str(name);
289 os.seekp(0, std::ios_base::end);
290 os << ": internal state #" << s % iNumDof << " of " << iNumDof
291 << " at point #" << s/iNumDof << " of " << GDI.iGetNum() << ", "
292 "block #" << s/(iNumDof*GDI.iGetNum()) << " of " << iNN;
293 desc[s] = os.str();
294 }
295
296 return;
297 }
298
299 desc.resize(1);
300 os << ": internal state #" << i % iNumDof << " of " << iNumDof
301 << " at point #" << i/iNumDof << " of " << GDI.iGetNum() << ", "
302 "block #" << i/(iNumDof*GDI.iGetNum()) << " of " << iNN;
303 desc[0] = os.str();
304 }
305
306 template <unsigned iNN>
307 std::ostream&
DescribeEq(std::ostream & out,const char * prefix,bool bInitial) const308 Aerodynamic2DElem<iNN>::DescribeEq(std::ostream& out,
309 const char *prefix, bool bInitial) const
310 {
311 ASSERT(bInitial == false);
312
313 integer iFirstIndex = iGetFirstIndex();
314 integer iNumDof = aerodata->iGetNumDof();
315
316 for (unsigned b = 0; b < iNN; b++) {
317 for (integer i = 0; i < GDI.iGetNum(); i++) {
318 out
319 << prefix << iFirstIndex + 1 << "->" << iFirstIndex + iNumDof
320 << ": dynamics equations at point #" << i << " of " << GDI.iGetNum() << ", "
321 "block #" << b << " of " << iNN << std::endl;
322
323 iFirstIndex += iNumDof;
324 }
325 }
326
327 return out;
328 }
329
330 template <unsigned iNN>
331 void
DescribeEq(std::vector<std::string> & desc,bool bInitial,int i) const332 Aerodynamic2DElem<iNN>::DescribeEq(std::vector<std::string>& desc,
333 bool bInitial, int i) const
334 {
335 if (i < -1) {
336 // error
337 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
338 }
339
340 ASSERT(bInitial == false);
341
342 std::ostringstream os;
343 os << elemnames[iNN - 1] << "(" << GetLabel() << ")";
344
345 integer iNumDof = aerodata->iGetNumDof();
346 if (i == -1) {
347 integer iTotDof = iNumDof*GDI.iGetNum();
348 desc.resize(iTotDof);
349
350 std::string name(os.str());
351
352 for (integer s = 0; s < iTotDof; s++) {
353 os.str(name);
354 os.seekp(0, std::ios_base::end);
355 os << ": internal equation #" << s % iNumDof << " of " << iNumDof
356 << " at point #" << s/iNumDof << " of " << GDI.iGetNum() << ", "
357 "block " << s/(iNumDof*GDI.iGetNum()) << " of " << iNN;
358 desc[s] = os.str();
359 }
360
361 return;
362 }
363
364 desc.resize(1);
365 os << ": internal equation #" << i % iNumDof << " of " << iNumDof
366 << " at point #" << i/iNumDof << " of " << GDI.iGetNum() << ", "
367 "block " << i/(iNumDof*GDI.iGetNum()) << " of " << iNN;
368 desc[0] = os.str();
369 }
370
371 template <unsigned iNN>
372 DofOrder::Order
GetDofType(unsigned int i) const373 Aerodynamic2DElem<iNN>::GetDofType(unsigned int i) const
374 {
375 ASSERT(aerodata->iGetNumDof() > 0);
376
377 return aerodata->GetDofType(i % aerodata->iGetNumDof());
378 }
379
380 template <unsigned iNN>
381 void
SetValue(DataManager * pDM,VectorHandler & X,VectorHandler & XP,SimulationEntity::Hints * h)382 Aerodynamic2DElem<iNN>::SetValue(DataManager *pDM,
383 VectorHandler& X, VectorHandler& XP,
384 SimulationEntity::Hints* h)
385 {
386 integer iNumDof = aerodata->iGetNumDof();
387 if (iNumDof) {
388 vx.Resize(6*iNN);
389 wx.Resize(6*iNN);
390 fq.Resize(iNumDof);
391 cq.Resize(iNumDof);
392
393 vx.Reset();
394 wx.Reset();
395 fq.Reset();
396 cq.Reset();
397 }
398 }
399
400 template <unsigned iNN>
401 void
AfterConvergence(const VectorHandler & X,const VectorHandler & XP)402 Aerodynamic2DElem<iNN>::AfterConvergence(
403 const VectorHandler& X,
404 const VectorHandler& XP)
405 {
406 /* Memoria in caso di forze instazionarie */
407 switch (aerodata->Unsteady()) {
408 case AeroData::STEADY:
409 break;
410
411 case AeroData::HARRIS:
412 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
413
414 case AeroData::BIELAWA:
415 for (unsigned i = 0; i < iNN*GDI.iGetNum(); i++) {
416 aerodata->Update(i);
417 }
418 break;
419
420 default:
421 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
422 }
423
424 for (unsigned i = 0; i < iNN*GDI.iGetNum(); i++) {
425 aerodata->AfterConvergence(i, X, XP);
426 }
427 }
428
429 /* Dimensioni del workspace */
430 template <unsigned iNN>
431 void
InitialWorkSpaceDim(integer * piNumRows,integer * piNumCols) const432 Aerodynamic2DElem<iNN>::InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const
433 {
434 *piNumRows = 6*iNN;
435 *piNumCols = 1;
436 }
437
438 /* assemblaggio jacobiano */
439 template <unsigned iNN>
440 VariableSubMatrixHandler&
InitialAssJac(VariableSubMatrixHandler & WorkMat,const VectorHandler &)441 Aerodynamic2DElem<iNN>::InitialAssJac(VariableSubMatrixHandler& WorkMat,
442 const VectorHandler& /* XCurr */)
443 {
444 DEBUGCOUTFNAME("Aerodynamic2DElem::InitialAssJac");
445 WorkMat.SetNullMatrix();
446 return WorkMat;
447 }
448
449 template <unsigned iNN>
450 void
OutputPrepare(OutputHandler & OH)451 Aerodynamic2DElem<iNN>::OutputPrepare(OutputHandler &OH)
452 {
453 if (bToBeOutput()) {
454 #ifdef USE_NETCDF
455 if (OH.UseNetCDF(OutputHandler::AERODYNAMIC)) {
456 ASSERT(OH.IsOpen(OutputHandler::NETCDF));
457
458 std::ostringstream os;
459 os << "elem.aerodynamic." << GetLabel();
460 (void)OH.CreateVar(os.str(), elemnames[iNN - 1]);
461
462 os << '.';
463 std::string name = os.str();
464
465 int totgp = iNN*GDI.iGetNum();
466 NetCDFOutputData.resize(totgp);
467
468 int j = 0;
469 for (std::vector<AeroNetCDFOutput>::iterator i = NetCDFOutputData.begin();
470 i != NetCDFOutputData.end(); ++i, ++j)
471 {
472 os.str("");
473 os << "Gauss point #" << j << "/" << totgp;
474 std::string gp(os.str());
475
476 unsigned uOutputFlags = (fToBeOutput() & OUTPUT_GP_ALL);
477
478 /* Add NetCDF (output) variables to the BinFile object
479 * and save the NcVar* pointer returned from add_var
480 * as handle for later write accesses.
481 * Define also variable attributes */
482 i->Var_X = 0;
483 if (uOutputFlags & AerodynamicOutput::OUTPUT_GP_X) {
484 os.str(name);
485 os.seekp(0, std::ios_base::end);
486 os << "X_" << j;
487 i->Var_X = OH.CreateVar<Vec3>(os.str(), "m",
488 gp + " global position vector (X, Y, Z)");
489 }
490
491 i->Var_Phi = 0;
492 if (uOutputFlags & AerodynamicOutput::OUTPUT_GP_R) {
493 os.str("");
494 os << "_" << j;
495 i->Var_Phi = OH.CreateRotationVar(name, os.str(), od,
496 gp + " global");
497 }
498
499 i->Var_V = 0;
500 if (uOutputFlags & AerodynamicOutput::OUTPUT_GP_V) {
501 os.str(name);
502 os.seekp(0, std::ios_base::end);
503 os << "V_" << j;
504 i->Var_V = OH.CreateVar<Vec3>(os.str(), "m/s",
505 gp + " global frame (V_X, V_Y, V_Z)");
506 }
507
508 i->Var_W = 0;
509 if (uOutputFlags & AerodynamicOutput::OUTPUT_GP_W) {
510 os.str(name);
511 os.seekp(0, std::ios_base::end);
512 os << "Omega_" << j;
513 i->Var_W = OH.CreateVar<Vec3>(os.str(), "radian/s",
514 gp + " global frame (Omega_X, Omega_Y, Omega_Z)");
515 }
516
517 i->Var_F = 0;
518 if (uOutputFlags & AerodynamicOutput::OUTPUT_GP_F) {
519 os.str(name);
520 os.seekp(0, std::ios_base::end);
521 os << "F_" << j;
522 i->Var_F = OH.CreateVar<Vec3>(os.str(), "N",
523 gp + " force in local frame (F_X, F_Y, F_Z)");
524 }
525
526 i->Var_M = 0;
527 if (uOutputFlags & AerodynamicOutput::OUTPUT_GP_M) {
528 os.str(name);
529 os.seekp(0, std::ios_base::end);
530 os << "M_" << j;
531 i->Var_M = OH.CreateVar<Vec3>(os.str(), "Nm",
532 gp + " force in local frame (M_X, M_Y, M_Z)");
533 }
534 }
535 }
536 #endif // USE_NETCDF
537 }
538 }
539
540 /* output; si assume che ogni tipo di elemento sappia, attraverso
541 * l'OutputHandler, dove scrivere il proprio output */
542 template <unsigned iNN>
543 void
Output_int(OutputHandler & OH) const544 Aerodynamic2DElem<iNN>::Output_int(OutputHandler &OH) const
545 {
546 #ifdef USE_NETCDF
547 if (OH.UseNetCDF(OutputHandler::AERODYNAMIC)) {
548 for (std::vector<AeroNetCDFOutput>::const_iterator i = NetCDFOutputData.begin();
549 i != NetCDFOutputData.end(); ++i)
550 {
551 if (i->Var_X) {
552 i->Var_X->put_rec(i->X.pGetVec(), OH.GetCurrentStep());
553 }
554
555 if (i->Var_Phi) {
556 Vec3 E;
557 switch (od) {
558 case EULER_123:
559 E = MatR2EulerAngles123(i->R)*dRaDegr;
560 break;
561
562 case EULER_313:
563 E = MatR2EulerAngles313(i->R)*dRaDegr;
564 break;
565
566 case EULER_321:
567 E = MatR2EulerAngles321(i->R)*dRaDegr;
568 break;
569
570 case ORIENTATION_VECTOR:
571 E = RotManip::VecRot(i->R);
572 break;
573
574 case ORIENTATION_MATRIX:
575 break;
576
577 default:
578 /* impossible */
579 break;
580 }
581
582 switch (od) {
583 case EULER_123:
584 case EULER_313:
585 case EULER_321:
586 case ORIENTATION_VECTOR:
587 i->Var_Phi->put_rec(E.pGetVec(), OH.GetCurrentStep());
588 break;
589
590 case ORIENTATION_MATRIX:
591 i->Var_Phi->put_rec(i->R.pGetMat(), OH.GetCurrentStep());
592 break;
593
594 default:
595 /* impossible */
596 break;
597 }
598 }
599
600 if (i->Var_V) {
601 i->Var_V->put_rec(i->V.pGetVec(), OH.GetCurrentStep());
602 }
603
604 if (i->Var_W) {
605 i->Var_W->put_rec(i->W.pGetVec(), OH.GetCurrentStep());
606 }
607
608 if (i->Var_F) {
609 i->Var_F->put_rec(i->F.pGetVec(), OH.GetCurrentStep());
610 }
611
612 if (i->Var_M) {
613 i->Var_M->put_rec(i->M.pGetVec(), OH.GetCurrentStep());
614 }
615 }
616 }
617 #endif /* USE_NETCDF */
618 }
619
620 // only send forces if:
621 // 1) an induced velocity model is defined
622 // 2) this element is not "passive" (i.e. contributes to induced velocity)
623 // 3) the induced velocity model does not require sectional forces
624 template <unsigned iNN>
625 void
AddForce_int(const StructNode * pN,const Vec3 & F,const Vec3 & M,const Vec3 & X) const626 Aerodynamic2DElem<iNN>::AddForce_int(const StructNode* pN,
627 const Vec3& F, const Vec3& M, const Vec3& X) const
628 {
629 if (pIndVel != 0 && !bPassiveInducedVelocity
630 && !pIndVel->bSectionalForces())
631 {
632 pIndVel->AddForce(this, pN, F, M, X);
633 }
634 }
635
636 // only send forces if:
637 // 1) an induced velocity model is defined
638 // 2) this element is not "passive" (i.e. contributes to induced velocity)
639 // 3) the induced velocity model requires sectional forces
640 template <unsigned iNN>
641 void
AddSectionalForce_int(unsigned uPnt,const Vec3 & F,const Vec3 & M,doublereal dW,const Vec3 & X,const Mat3x3 & R,const Vec3 & V,const Vec3 & W) const642 Aerodynamic2DElem<iNN>::AddSectionalForce_int(unsigned uPnt,
643 const Vec3& F, const Vec3& M, doublereal dW,
644 const Vec3& X, const Mat3x3& R,
645 const Vec3& V, const Vec3& W) const
646 {
647 if (pIndVel != 0 && !bPassiveInducedVelocity
648 && pIndVel->bSectionalForces())
649 {
650 pIndVel->AddSectionalForce(GetElemType(), this, uPnt,
651 F, M, dW, X, R, V, W);
652 }
653 }
654
655 /* Aerodynamic2DElem - end */
656
657
658 /* AerodynamicBody - begin */
659
AerodynamicBody(unsigned int uLabel,const DofOwner * pDO,const StructNode * pN,InducedVelocity * pR,bool bPassive,const Vec3 & fTmp,doublereal dS,const Mat3x3 & RaTmp,const Shape * pC,const Shape * pF,const Shape * pV,const Shape * pT,const Shape * pTL,integer iN,AeroData * a,const DriveCaller * pDC,bool bUseJacobian,OrientationDescription ood,flag fOut)660 AerodynamicBody::AerodynamicBody(unsigned int uLabel,
661 const DofOwner *pDO,
662 const StructNode* pN, InducedVelocity* pR, bool bPassive,
663 const Vec3& fTmp, doublereal dS,
664 const Mat3x3& RaTmp,
665 const Shape* pC, const Shape* pF,
666 const Shape* pV, const Shape* pT,
667 const Shape* pTL,
668 integer iN, AeroData* a,
669 const DriveCaller* pDC,
670 bool bUseJacobian,
671 OrientationDescription ood,
672 flag fOut)
673 : Elem(uLabel, fOut),
674 Aerodynamic2DElem<1>(uLabel, pDO, pR, bPassive, pC, pF, pV, pT, pTL, iN,
675 a, pDC, bUseJacobian, ood, fOut),
676 pNode(pN),
677 f(fTmp),
678 dHalfSpan(dS/2.),
679 Ra(RaTmp),
680 Ra3(RaTmp.GetVec(3)),
681 F(Zero3),
682 M(Zero3)
683 {
684 DEBUGCOUTFNAME("AerodynamicBody::AerodynamicBody");
685
686 ASSERT(pNode != 0);
687 ASSERT(pNode->GetNodeType() == Node::STRUCTURAL);
688
689 }
690
~AerodynamicBody(void)691 AerodynamicBody::~AerodynamicBody(void)
692 {
693 DEBUGCOUTFNAME("AerodynamicBody::~AerodynamicBody");
694 }
695
696 /* Scrive il contributo dell'elemento al file di restart */
697 std::ostream&
Restart(std::ostream & out) const698 AerodynamicBody::Restart(std::ostream& out) const
699 {
700 DEBUGCOUTFNAME("AerodynamicBody::Restart");
701
702 out << " aerodynamic body: " << GetLabel() << ", "
703 << pNode->GetLabel();
704 if (pIndVel != 0) {
705 out << ", rotor, " << pIndVel->GetLabel();
706 }
707 out << ", reference, node, ", f.Write(out, ", ")
708 << ", " << dHalfSpan*2.
709 << ", reference, node, 1, ", (Ra.GetVec(1)).Write(out, ", ")
710 << ", 2, ", (Ra.GetVec(2)).Write(out, ", ")
711 << ", ";
712 Chord.pGetShape()->Restart(out) << ", ";
713 ForcePoint.pGetShape()->Restart(out) << ", ";
714 VelocityPoint.pGetShape()->Restart(out) << ", ";
715 Twist.pGetShape()->Restart(out) << ", "
716 << ", " << GDI.iGetNum() << ", control, ";
717 pGetDriveCaller()->Restart(out) << ", ";
718 aerodata->Restart(out);
719 return out << ";" << std::endl;
720 }
721
722 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)723 AerodynamicBody::AssJac(VariableSubMatrixHandler& WorkMat,
724 doublereal dCoef,
725 const VectorHandler& XCurr,
726 const VectorHandler& XPrimeCurr)
727 {
728 DEBUGCOUT("Entering AerodynamicBody::AssJac()" << std::endl);
729
730 if (!bJacobian) {
731 WorkMat.SetNullMatrix();
732 return WorkMat;
733 }
734
735 FullSubMatrixHandler& WM = WorkMat.SetFull();
736
737 /* Ridimensiona la sottomatrice in base alle esigenze */
738 integer iNumRows = 0;
739 integer iNumCols = 0;
740 WorkSpaceDim(&iNumRows, &iNumCols);
741 WM.ResizeReset(iNumRows, iNumCols);
742
743 /* Recupera gli indici delle varie incognite */
744 integer iNodeFirstMomIndex = pNode->iGetFirstMomentumIndex();
745 integer iNodeFirstPosIndex = pNode->iGetFirstPositionIndex();
746
747 /* Setta gli indici delle equazioni */
748 for (int iCnt = 1; iCnt <= 6; iCnt++) {
749 WM.PutRowIndex(iCnt, iNodeFirstMomIndex + iCnt);
750 WM.PutColIndex(iCnt, iNodeFirstPosIndex + iCnt);
751 }
752
753 /* Equations start here... */
754 doublereal dW[6];
755
756 /* Dati del nodo */
757 const Vec3& Xn(pNode->GetXCurr());
758 const Mat3x3& Rn(pNode->GetRCurr());
759 const Vec3& Vn(pNode->GetVCurr());
760 const Vec3& Wn(pNode->GetWCurr());
761
762 /*
763 * Matrice di trasformazione dal sistema globale a quello aerodinamico
764 */
765 Mat3x3 RR(Rn*Ra);
766
767 /*
768 * Se l'elemento e' collegato ad un rotore,
769 * si fa dare la velocita' di rotazione
770 */
771 doublereal dOmega = 0.;
772 if (pIndVel != 0) {
773 Rotor *pRotor = dynamic_cast<Rotor *>(pIndVel);
774 if (pRotor) {
775 dOmega = pRotor->dGetOmega();
776 }
777 }
778
779 /*
780 * Dati "permanenti" (uso la posizione del nodo perche'
781 * non dovrebbero cambiare "molto")
782 */
783 doublereal rho, c, p, T;
784 GetAirProps(Xn, rho, c, p, T); /* p, T no used yet */
785 aerodata->SetAirData(rho, c);
786
787 ResetIterator();
788
789 integer iNumDof = aerodata->iGetNumDof();
790 integer iFirstEq = -1;
791 integer iFirstSubEq = -1;
792 if (iNumDof > 0) {
793 iFirstEq = iGetFirstIndex();
794 iFirstSubEq = 6;
795
796 integer iOffset = iFirstEq - 6;
797
798 for (int iCnt = 6 + 1; iCnt <= iNumRows; iCnt++) {
799 WM.PutRowIndex(iCnt, iOffset + iCnt);
800 WM.PutColIndex(iCnt, iOffset + iCnt);
801 }
802 }
803
804 /* Ciclo sui punti di Gauss */
805 PntWght PW = GDI.GetFirst();
806 int iPnt = 0;
807 do {
808 doublereal dCsi = PW.dGetPnt();
809 Vec3 Xr(Rn*(f + Ra3*(dHalfSpan*dCsi)));
810 Vec3 Xnr = Xn + Xr;
811 Vec3 Vr(Vn + Wn.Cross(Xr));
812
813 /* Contributo di velocita' del vento */
814 Vec3 VTmp(Zero3);
815 if (fGetAirVelocity(VTmp, Xnr)) {
816 Vr -= VTmp;
817 }
818
819 /*
820 * Se l'elemento e' collegato ad un rotore,
821 * aggiunge alla velocita' la velocita' indotta
822 */
823 if (pIndVel != 0) {
824 Vr += pIndVel->GetInducedVelocity(GetElemType(),
825 GetLabel(), iPnt, Xnr);
826 }
827
828 /* Copia i dati nel vettore di lavoro dVAM */
829 doublereal dTw = Twist.dGet(dCsi) + dGet();
830 aerodata->SetSectionData(dCsi,
831 Chord.dGet(dCsi),
832 ForcePoint.dGet(dCsi),
833 VelocityPoint.dGet(dCsi),
834 dTw,
835 dOmega);
836
837 /*
838 * Lo svergolamento non viene piu' trattato in aerod2_; quindi
839 * lo uso per correggere la matrice di rotazione
840 * dal sistema aerodinamico a quello globale
841 */
842 Mat3x3 RRloc;
843 if (dTw != 0.) {
844 doublereal dCosT = cos(dTw);
845 doublereal dSinT = sin(dTw);
846 /* Assumo lo svergolamento positivo a cabrare */
847 Mat3x3 RTw( dCosT, dSinT, 0.,
848 -dSinT, dCosT, 0.,
849 0., 0., 1.);
850 RRloc = RR*RTw;
851
852 } else {
853 RRloc = RR;
854 }
855
856 /*
857 * Ruota velocita' e velocita' angolare nel sistema
858 * aerodinamico e li copia nel vettore di lavoro dW
859 */
860 VTmp = RRloc.MulTV(Vr);
861 VTmp.PutTo(dW);
862
863 Vec3 WTmp = RRloc.MulTV(Wn);
864 WTmp.PutTo(&dW[3]);
865
866 /* Funzione di calcolo delle forze aerodinamiche */
867 doublereal Fa0[6];
868 Mat6x6 JFa;
869
870 /* Jacobian Assembly... */
871 doublereal cc = dHalfSpan*PW.dGetWght();
872
873 if (iNumDof) {
874 // prepare (v/dot{x} + dCoef*v/x) and so
875 Mat3x3 RRlocT(RRloc.Transpose());
876
877 vx.PutMat3x3(1, RRlocT);
878 vx.PutMat3x3(4, RRloc.MulTM(Mat3x3(MatCross, (Vr + Xr.Cross(Wn))*dCoef - Xr)));
879 wx.PutMat3x3(4, RRlocT);
880
881 // equations from iFirstEq on are dealt with by aerodata
882 aerodata->AssJac(WM, dCoef, XCurr, XPrimeCurr,
883 iFirstEq, iFirstSubEq,
884 vx, wx, fq, cq, iPnt, dW, Fa0, JFa, OUTA[iPnt]);
885
886 // deal with (f/dot{q} + dCoef*f/q) and so
887 integer iOffset = 6 + iPnt*iNumDof;
888 for (integer iCol = 1; iCol <= iNumDof; iCol++) {
889 Vec3 fqTmp((RRloc*fq.GetVec(iCol))*cc);
890 Vec3 cqTmp(Xr.Cross(fqTmp) + (RRloc*cq.GetVec(iCol))*cc);
891
892 WM.Sub(1, iOffset + iCol, fqTmp);
893 WM.Sub(4, iOffset + iCol, cqTmp);
894 }
895
896 // first equation
897 iFirstEq += iNumDof;
898 iFirstSubEq += iNumDof;
899
900 } else {
901 aerodata->GetForcesJac(iPnt, dW, Fa0, JFa, OUTA[iPnt]);
902 }
903
904 // rotate force, couple and Jacobian matrix in absolute frame
905 Mat6x6 JFaR = MultRMRt(JFa, RRloc, cc);
906
907 Vec3 fTmp(RRloc*(Vec3(&Fa0[0])*dCoef));
908 Vec3 cTmp(RRloc*(Vec3(&Fa0[3])*dCoef) + Xr.Cross(fTmp));
909
910 // compute submatrices (see tecman.pdf)
911 Mat3x3 Mat21(Xr.Cross(JFaR.GetMat11()) + JFaR.GetMat21());
912
913 Mat3x3 Mat12(JFaR.GetMat12() - JFaR.GetMat11()*Mat3x3(MatCross, Xr));
914
915 Mat3x3 Mat22(Xr.Cross(JFaR.GetMat12()) + JFaR.GetMat22() - Mat21*Mat3x3(MatCross, Xr));
916
917 Mat3x3 MatV(MatCross, (Vr + Xr.Cross(Wn))*dCoef);
918
919 Mat12 += JFaR.GetMat11()*MatV - Mat3x3(MatCross, fTmp);
920
921 Mat22 += Mat21*MatV - Mat3x3(MatCross, cTmp);
922
923 // add (actually, sub) contributions, scaled by weight
924 WM.Sub(1, 1, JFaR.GetMat11());
925 WM.Sub(1, 4, Mat12);
926 WM.Sub(4, 1, Mat21);
927 WM.Sub(4, 4, Mat22);
928
929 iPnt++;
930
931 } while (GDI.fGetNext(PW));
932
933 return WorkMat;
934 }
935
936
937 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)938 AerodynamicBody::AssRes(SubVectorHandler& WorkVec,
939 doublereal dCoef,
940 const VectorHandler& XCurr,
941 const VectorHandler& XPrimeCurr)
942 {
943 DEBUGCOUTFNAME("AerodynamicBody::AssRes");
944
945 integer iNumRows;
946 integer iNumCols;
947 WorkSpaceDim(&iNumRows, &iNumCols);
948 WorkVec.ResizeReset(iNumRows);
949
950 integer iFirstIndex = pNode->iGetFirstMomentumIndex();
951 for (int iCnt = 1; iCnt <= 6; iCnt++) {
952 WorkVec.PutRowIndex(iCnt, iFirstIndex + iCnt);
953 }
954
955 AssVec(WorkVec, dCoef, XCurr, XPrimeCurr);
956
957 return WorkVec;
958 }
959
960
961 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)962 AerodynamicBody::InitialAssRes(SubVectorHandler& WorkVec,
963 const VectorHandler& XCurr)
964 {
965 DEBUGCOUTFNAME("AerodynamicBody::InitialAssRes");
966
967 if (aerodata->iGetNumDof() > 0) {
968 WorkVec.ResizeReset(0);
969 return WorkVec;
970 }
971
972 integer iNumRows;
973 integer iNumCols;
974 WorkSpaceDim(&iNumRows, &iNumCols);
975 WorkVec.ResizeReset(iNumRows);
976
977 integer iFirstIndex = pNode->iGetFirstPositionIndex();
978 for (int iCnt = 1; iCnt <= 6; iCnt++) {
979 WorkVec.PutRowIndex(iCnt, iFirstIndex + iCnt);
980 }
981
982 AssVec(WorkVec, 1., XCurr, XCurr);
983
984 return WorkVec;
985 }
986
987
988 /* assemblaggio residuo */
989 void
AssVec(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)990 AerodynamicBody::AssVec(SubVectorHandler& WorkVec,
991 doublereal dCoef,
992 const VectorHandler& XCurr,
993 const VectorHandler& XPrimeCurr)
994 {
995 DEBUGCOUTFNAME("AerodynamicBody::AssVec");
996
997 doublereal dTng[6];
998 doublereal dW[6];
999
1000 /* Dati del nodo */
1001 const Vec3& Xn(pNode->GetXCurr());
1002 const Mat3x3& Rn(pNode->GetRCurr());
1003 const Vec3& Vn(pNode->GetVCurr());
1004 const Vec3& Wn(pNode->GetWCurr());
1005
1006 /*
1007 * Matrice di trasformazione dal sistema globale a quello aerodinamico
1008 */
1009 Mat3x3 RR(Rn*Ra);
1010
1011 /*
1012 * Se l'elemento e' collegato ad un rotore,
1013 * si fa dare la velocita' di rotazione
1014 */
1015 doublereal dOmega = 0.;
1016 if (pIndVel != 0) {
1017 Rotor *pRotor = dynamic_cast<Rotor *>(pIndVel);
1018 if (pRotor) {
1019 dOmega = pRotor->dGetOmega();
1020 }
1021 }
1022
1023 /* Resetta i dati */
1024 F.Reset();
1025 M.Reset();
1026
1027 /*
1028 * Dati "permanenti" (uso la posizione del nodo perche'
1029 * non dovrebbero cambiare "molto")
1030 */
1031 doublereal rho, c, p, T;
1032 GetAirProps(Xn, rho, c, p, T); /* p, T no used yet */
1033 aerodata->SetAirData(rho, c);
1034
1035 ResetIterator();
1036
1037 integer iNumDof = aerodata->iGetNumDof();
1038 integer iFirstEq = -1;
1039 integer iFirstSubEq = -1;
1040 if (iNumDof > 0) {
1041 iFirstEq = iGetFirstIndex();
1042 iFirstSubEq = 6;
1043
1044 integer iOffset = iFirstEq - 6;
1045 integer iNumRows = WorkVec.iGetSize();
1046 for (int iCnt = 6 + 1; iCnt <= iNumRows; iCnt++) {
1047 WorkVec.PutRowIndex(iCnt, iOffset + iCnt);
1048 }
1049 }
1050
1051 /* Ciclo sui punti di Gauss */
1052 PntWght PW = GDI.GetFirst();
1053 int iPnt = 0;
1054 do {
1055 doublereal dCsi = PW.dGetPnt();
1056 Vec3 Xr(Rn*(f + Ra3*(dHalfSpan*dCsi)));
1057 Vec3 Xnr = Xn + Xr;
1058 Vec3 Vr(Vn + Wn.Cross(Xr));
1059
1060 /* Contributo di velocita' del vento */
1061 Vec3 VTmp(Zero3);
1062 if (fGetAirVelocity(VTmp, Xnr)) {
1063 Vr -= VTmp;
1064 }
1065
1066 /*
1067 * Se l'elemento e' collegato ad un rotore,
1068 * aggiunge alla velocita' la velocita' indotta
1069 */
1070 if (pIndVel != 0) {
1071 Vr += pIndVel->GetInducedVelocity(GetElemType(),
1072 GetLabel(), iPnt, Xnr);
1073 }
1074
1075 /* Copia i dati nel vettore di lavoro dVAM */
1076 doublereal dTw = Twist.dGet(dCsi) + dGet();
1077 aerodata->SetSectionData(dCsi,
1078 Chord.dGet(dCsi),
1079 ForcePoint.dGet(dCsi),
1080 VelocityPoint.dGet(dCsi),
1081 dTw,
1082 dOmega);
1083
1084 /*
1085 * Lo svergolamento non viene piu' trattato in aerod2_; quindi
1086 * lo uso per correggere la matrice di rotazione
1087 * dal sistema aerodinamico a quello globale
1088 */
1089 Mat3x3 RRloc;
1090 if (dTw != 0.) {
1091 doublereal dCosT = cos(dTw);
1092 doublereal dSinT = sin(dTw);
1093 /* Assumo lo svergolamento positivo a cabrare */
1094 Mat3x3 RTw( dCosT, dSinT, 0.,
1095 -dSinT, dCosT, 0.,
1096 0., 0., 1.);
1097 RRloc = RR*RTw;
1098
1099 } else {
1100 RRloc = RR;
1101 }
1102
1103 /*
1104 * Ruota velocita' e velocita' angolare nel sistema
1105 * aerodinamico e li copia nel vettore di lavoro dW
1106 */
1107 VTmp = RRloc.MulTV(Vr);
1108 VTmp.PutTo(dW);
1109
1110 Vec3 WTmp = RRloc.MulTV(Wn);
1111 WTmp.PutTo(&dW[3]);
1112
1113 /* Funzione di calcolo delle forze aerodinamiche */
1114 if (iNumDof) {
1115 aerodata->AssRes(WorkVec, dCoef, XCurr, XPrimeCurr,
1116 iFirstEq, iFirstSubEq, iPnt, dW, dTng, OUTA[iPnt]);
1117
1118 // first equation
1119 iFirstEq += iNumDof;
1120 iFirstSubEq += iNumDof;
1121
1122 } else {
1123 aerodata->GetForces(iPnt, dW, dTng, OUTA[iPnt]);
1124 }
1125
1126 /* Dimensionalizza le forze */
1127 doublereal dWght = dHalfSpan*PW.dGetWght();
1128 dTng[1] *= TipLoss.dGet(dCsi);
1129 Vec3 FTmp(RRloc*(Vec3(&dTng[0])));
1130 Vec3 MTmp(RRloc*(Vec3(&dTng[3])));
1131
1132 // Se e' definito il rotore, aggiungere il contributo alla trazione
1133 AddSectionalForce_int(iPnt, FTmp, MTmp, dWght, Xnr, RRloc, Vr, Wn);
1134
1135 FTmp *= dWght;
1136 MTmp *= dWght;
1137
1138 F += FTmp;
1139 M += MTmp;
1140 M += Xr.Cross(FTmp);
1141
1142 // specific for Gauss points force output
1143 if (bToBeOutput()) {
1144 SetData(VTmp, dTng, Xr, RRloc, Vr, Wn, FTmp, MTmp);
1145 }
1146
1147 iPnt++;
1148
1149 } while (GDI.fGetNext(PW));
1150
1151 // Se e' definito il rotore, aggiungere il contributo alla trazione
1152 AddForce_int(pNode, F, M, Xn);
1153
1154 /* Sommare il termine al residuo */
1155 WorkVec.Add(1, F);
1156 WorkVec.Add(4, M);
1157 }
1158
1159 /*
1160 * output; si assume che ogni tipo di elemento sappia, attraverso
1161 * l'OutputHandler, dove scrivere il proprio output
1162 */
1163 void
Output(OutputHandler & OH) const1164 AerodynamicBody::Output(OutputHandler& OH) const
1165 {
1166 /* Output delle forze aerodinamiche F, M su apposito file */
1167 if (bToBeOutput()) {
1168 Aerodynamic2DElem<1>::Output_int(OH);
1169
1170 if (OH.UseText(OutputHandler::AERODYNAMIC)) {
1171 std::ostream& out = OH.Aerodynamic()
1172 << std::setw(8) << GetLabel();
1173
1174 switch (GetOutput()) {
1175 case AEROD_OUT_NODE:
1176 out << " " << std::setw(8) << pNode->GetLabel()
1177 << " ", F.Write(out, " ") << " ", M.Write(out, " ");
1178 break;
1179
1180 case AEROD_OUT_PGAUSS:
1181 ASSERT(!OutputData.empty());
1182 for (std::vector<Aero_output>::const_iterator i = OutputData.begin();
1183 i != OutputData.end(); ++i)
1184 {
1185 out << " " << i->alpha
1186 << " " << i->f;
1187 }
1188 break;
1189
1190 case AEROD_OUT_STD:
1191 for (int i = 0; i < GDI.iGetNum(); i++) {
1192 out
1193 << " " << OUTA[i].alpha
1194 << " " << OUTA[i].gamma
1195 << " " << OUTA[i].mach
1196 << " " << OUTA[i].cl
1197 << " " << OUTA[i].cd
1198 << " " << OUTA[i].cm
1199 << " " << OUTA[i].alf1
1200 << " " << OUTA[i].alf2;
1201 }
1202 break;
1203
1204 default:
1205 ASSERT(0);
1206 break;
1207 }
1208
1209 out << std::endl;
1210 } }
1211 }
1212
1213 /* AerodynamicBody - end */
1214
1215 static bool
ReadInducedVelocity(DataManager * pDM,MBDynParser & HP,unsigned uLabel,const char * sElemType,InducedVelocity * & pIndVel,bool & bPassive)1216 ReadInducedVelocity(DataManager *pDM, MBDynParser& HP,
1217 unsigned uLabel, const char *sElemType,
1218 InducedVelocity*& pIndVel, bool &bPassive)
1219 {
1220 bool bReadIV(false);
1221 bool bReadUDIV(false);
1222 if (HP.IsKeyWord("rotor")) {
1223 silent_cerr(sElemType << "(" << uLabel << "): "
1224 "\"rotor\" keyword is deprecated; "
1225 "use \"induced velocity\" instead "
1226 "at line " << HP.GetLineData()
1227 << std::endl);
1228
1229 bReadIV = true;
1230
1231 } else if (HP.IsKeyWord("induced" "velocity")) {
1232 bReadIV = true;
1233
1234 } else if (HP.IsKeyWord("user" "defined" "induced" "velocity")) {
1235 bReadIV = true;
1236 bReadUDIV = true;
1237 }
1238
1239 if (bReadIV) {
1240 unsigned int uIV = (unsigned int)HP.GetInt();
1241 DEBUGLCOUT(MYDEBUG_INPUT,
1242 "Linked to InducedVelocity(" << uIV << ")" << std::endl);
1243
1244 bPassive = false;
1245 if (HP.IsKeyWord("passive")) {
1246 bPassive = true;
1247 }
1248
1249 /*
1250 * verifica di esistenza del rotore
1251 * NOTA: ovviamente il rotore deve essere definito
1252 * prima dell'elemento aerodinamico
1253 */
1254 Elem* p;
1255
1256 if (bReadUDIV) {
1257 p = pDM->pFindElem(Elem::LOADABLE, uIV);
1258 if (p == 0) {
1259 silent_cerr(sElemType << "(" << uLabel << "): "
1260 "user-defined InducedVelocity(" << uIV << ") not defined "
1261 "at line " << HP.GetLineData()
1262 << std::endl);
1263 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1264 }
1265
1266 } else {
1267 p = pDM->pFindElem(Elem::INDUCEDVELOCITY, uIV);
1268 if (p == 0) {
1269 // try a user-defined one?
1270 p = pDM->pFindElem(Elem::LOADABLE, uIV);
1271 if (p == 0 || !dynamic_cast<InducedVelocity *>(p)) {
1272 silent_cerr(sElemType << "(" << uLabel << "): "
1273 "InducedVelocity(" << uIV << ") not defined "
1274 "at line " << HP.GetLineData()
1275 << std::endl);
1276
1277 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1278 }
1279
1280 silent_cerr(sElemType << "(" << uLabel << "): "
1281 "InducedVelocity(" << uIV << ") not defined; using user-defined InducedVelocity(" << uIV << ") "
1282 "at line " << HP.GetLineData()
1283 << std::endl);
1284 }
1285 }
1286
1287 pIndVel = dynamic_cast<InducedVelocity *>(p);
1288 ASSERT(pIndVel != 0);
1289 }
1290
1291 return (pIndVel != 0);
1292 }
1293
1294 void
ReadAerodynamicCustomOutput(DataManager * pDM,MBDynParser & HP,unsigned int uLabel,unsigned & uFlags,OrientationDescription & od)1295 ReadAerodynamicCustomOutput(DataManager* pDM, MBDynParser& HP, unsigned int uLabel,
1296 unsigned& uFlags, OrientationDescription& od)
1297 {
1298 uFlags = AerodynamicOutput::OUTPUT_NONE;
1299
1300 while (HP.IsArg()) {
1301 unsigned uFlag;
1302
1303 if (HP.IsKeyWord("position")) {
1304 uFlag = AerodynamicOutput::OUTPUT_GP_X;
1305
1306 } else if (HP.IsKeyWord("orientation")) {
1307 uFlag = AerodynamicOutput::OUTPUT_GP_R;
1308
1309 } else if (HP.IsKeyWord("velocity")) {
1310 uFlag = AerodynamicOutput::OUTPUT_GP_V;
1311
1312 } else if (HP.IsKeyWord("angular" "velocity")) {
1313 uFlag = AerodynamicOutput::OUTPUT_GP_W;
1314
1315 } else if (HP.IsKeyWord("configuration")) {
1316 uFlag = AerodynamicOutput::OUTPUT_GP_CONFIGURATION;
1317
1318 } else if (HP.IsKeyWord("force")) {
1319 uFlag = AerodynamicOutput::OUTPUT_GP_F;
1320
1321 } else if (HP.IsKeyWord("moment")) {
1322 uFlag = AerodynamicOutput::OUTPUT_GP_M;
1323
1324 } else if (HP.IsKeyWord("forces")) {
1325 uFlag = AerodynamicOutput::OUTPUT_GP_FORCES;
1326
1327 } else if (HP.IsKeyWord("all")) {
1328 uFlag = AerodynamicOutput::OUTPUT_GP_ALL;
1329
1330 } else {
1331 break;
1332 }
1333
1334 if (uFlags & uFlag) {
1335 silent_cerr("AerodynamicElement(" << uLabel << "): "
1336 "duplicate custom output "
1337 "at line " << HP.GetLineData()
1338 << std::endl);
1339 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1340 }
1341
1342 if (uFlag & AerodynamicOutput::OUTPUT_GP_R) {
1343 od = ReadOptionalOrientationDescription(pDM, HP);
1344 }
1345
1346 uFlags |= uFlag;
1347 }
1348 }
1349
1350 void
ReadOptionalAerodynamicCustomOutput(DataManager * pDM,MBDynParser & HP,unsigned int uLabel,unsigned & uFlags,OrientationDescription & od)1351 ReadOptionalAerodynamicCustomOutput(DataManager* pDM, MBDynParser& HP, unsigned int uLabel,
1352 unsigned& uFlags, OrientationDescription& od)
1353 {
1354 pDM->GetOutput(Elem::AERODYNAMIC, uFlags, od);
1355 if (HP.IsKeyWord("custom" "output")) {
1356 ReadAerodynamicCustomOutput(pDM, HP, uLabel, uFlags, od);
1357 }
1358 }
1359
1360 Elem *
ReadAerodynamicBody(DataManager * pDM,MBDynParser & HP,const DofOwner * pDO,unsigned int uLabel)1361 ReadAerodynamicBody(DataManager* pDM,
1362 MBDynParser& HP,
1363 const DofOwner *pDO,
1364 unsigned int uLabel)
1365 {
1366 DEBUGCOUTFNAME("ReadAerodynamicBody");
1367
1368 /* Nodo */
1369 const StructNode* pNode = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
1370
1371 InducedVelocity* pIndVel = 0;
1372 bool bPassive(false);
1373 (void)ReadInducedVelocity(pDM, HP, uLabel, "AerodynamicBody",
1374 pIndVel, bPassive);
1375
1376 ReferenceFrame RF(pNode);
1377 Vec3 f(HP.GetPosRel(RF));
1378
1379 DEBUGLCOUT(MYDEBUG_INPUT, "Offset: " << f << std::endl);
1380
1381 Mat3x3 Ra(HP.GetRotRel(RF));
1382
1383 doublereal dSpan = HP.GetReal();
1384 DEBUGLCOUT(MYDEBUG_INPUT, "Span: " << dSpan << std::endl);
1385
1386 Shape* pChord = 0;
1387 Shape* pForce = 0;
1388 Shape* pVelocity = 0;
1389 Shape* pTwist = 0;
1390 Shape* pTipLoss = 0;
1391
1392 integer iNumber = 0;
1393 DriveCaller* pDC = 0;
1394 AeroData* aerodata = 0;
1395
1396 ReadAeroData(pDM, HP, 1,
1397 &pChord, &pForce, &pVelocity, &pTwist, &pTipLoss,
1398 &iNumber, &pDC, &aerodata);
1399
1400 bool bUseJacobian(false);
1401 if (HP.IsKeyWord("jacobian")) {
1402 bUseJacobian = HP.GetYesNoOrBool(bDefaultUseJacobian);
1403 }
1404
1405 if (aerodata->iGetNumDof() > 0 && !bUseJacobian) {
1406 silent_cerr("AerodynamicBody(" << uLabel << "): "
1407 "aerodynamic model needs \"jacobian, yes\" at line " << HP.GetLineData()
1408 << std::endl);
1409 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1410 }
1411
1412 OrientationDescription od = UNKNOWN_ORIENTATION_DESCRIPTION;
1413 unsigned uFlags = AerodynamicOutput::OUTPUT_NONE;
1414 ReadOptionalAerodynamicCustomOutput(pDM, HP, uLabel, uFlags, od);
1415
1416 flag fOut = pDM->fReadOutput(HP, Elem::AERODYNAMIC);
1417 if (HP.IsArg()) {
1418 if (HP.IsKeyWord("std")) {
1419 fOut |= AerodynamicOutput::AEROD_OUT_STD;
1420 } else if (HP.IsKeyWord("gauss")) {
1421 fOut |= AerodynamicOutput::AEROD_OUT_PGAUSS;
1422 } else if (HP.IsKeyWord("node")) {
1423 fOut |= AerodynamicOutput::AEROD_OUT_NODE;
1424 } else {
1425 silent_cerr("AerodynamicBody(" << uLabel << "): "
1426 "unknown output mode at line " << HP.GetLineData()
1427 << std::endl);
1428 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1429 }
1430
1431 } else if (fOut) {
1432 fOut |= AerodynamicOutput::AEROD_OUT_STD;
1433 }
1434 fOut |= uFlags;
1435
1436 Elem* pEl = 0;
1437 SAFENEWWITHCONSTRUCTOR(pEl,
1438 AerodynamicBody,
1439 AerodynamicBody(uLabel, pDO, pNode, pIndVel, bPassive,
1440 f, dSpan, Ra,
1441 pChord, pForce, pVelocity, pTwist, pTipLoss,
1442 iNumber, aerodata, pDC, bUseJacobian, od, fOut));
1443
1444 /* Se non c'e' il punto e virgola finale */
1445 if (HP.IsArg()) {
1446 silent_cerr("semicolon expected at line "
1447 << HP.GetLineData() << std::endl);
1448 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1449 }
1450
1451 Vec3 Ra3 = Ra.GetVec(3);
1452 doublereal dCm1 = pChord->dGet(-1.);
1453 doublereal dPm1 = pForce->dGet(-1.);
1454 doublereal dTm1 = pTwist->dGet(-1.);
1455 Vec3 Ram1 = Ra*Vec3(std::cos(dTm1), std::sin(dTm1), 0.);
1456
1457 doublereal dCp1 = pChord->dGet(1.);
1458 doublereal dPp1 = pForce->dGet(1.);
1459 doublereal dTp1 = pTwist->dGet(1.);
1460 Vec3 Rap1 = Ra*Vec3(std::cos(dTp1), std::sin(dTp1), 0.);
1461
1462 std::ostream& out = pDM->GetLogFile();
1463 out << "aero0: " << uLabel
1464 << " " << pNode->GetLabel()
1465 << " ", (f - Ra3*(dSpan/2.) + Ram1*(dPm1 - dCm1*3./4.)).Write(out, " ")
1466 << " ", (f - Ra3*(dSpan/2.) + Ram1*(dPm1 + dCm1/4.)).Write(out, " ")
1467 << " ", (f + Ra3*(dSpan/2.) + Rap1*(dPp1 - dCp1*3./4.)).Write(out, " ")
1468 << " ", (f + Ra3*(dSpan/2.) + Rap1*(dPp1 + dCp1/4.)).Write(out, " ")
1469 << std::endl;
1470
1471 return pEl;
1472 } /* End of ReadAerodynamicBody() */
1473
1474
1475 /* AerodynamicBeam - begin */
1476
AerodynamicBeam(unsigned int uLabel,const DofOwner * pDO,const Beam * pB,InducedVelocity * pR,bool bPassive,const Vec3 & fTmp1,const Vec3 & fTmp2,const Vec3 & fTmp3,const Mat3x3 & Ra1Tmp,const Mat3x3 & Ra2Tmp,const Mat3x3 & Ra3Tmp,const Shape * pC,const Shape * pF,const Shape * pV,const Shape * pT,const Shape * pTL,integer iN,AeroData * a,const DriveCaller * pDC,bool bUseJacobian,OrientationDescription ood,flag fOut)1477 AerodynamicBeam::AerodynamicBeam(unsigned int uLabel,
1478 const DofOwner *pDO,
1479 const Beam* pB, InducedVelocity* pR, bool bPassive,
1480 const Vec3& fTmp1,
1481 const Vec3& fTmp2,
1482 const Vec3& fTmp3,
1483 const Mat3x3& Ra1Tmp,
1484 const Mat3x3& Ra2Tmp,
1485 const Mat3x3& Ra3Tmp,
1486 const Shape* pC, const Shape* pF,
1487 const Shape* pV, const Shape* pT,
1488 const Shape* pTL,
1489 integer iN, AeroData* a,
1490 const DriveCaller* pDC,
1491 bool bUseJacobian,
1492 OrientationDescription ood,
1493 flag fOut)
1494 : Elem(uLabel, fOut),
1495 Aerodynamic2DElem<3>(uLabel, pDO, pR, bPassive,
1496 pC, pF, pV, pT, pTL, iN, a, pDC, bUseJacobian, ood, fOut),
1497 pBeam(pB),
1498 f1(fTmp1),
1499 f2(fTmp2),
1500 f3(fTmp3),
1501 Ra1(Ra1Tmp),
1502 Ra2(Ra2Tmp),
1503 Ra3(Ra3Tmp),
1504 Ra1_3(Ra1Tmp.GetVec(3)),
1505 Ra2_3(Ra2Tmp.GetVec(3)),
1506 Ra3_3(Ra3Tmp.GetVec(3))
1507 {
1508 DEBUGCOUTFNAME("AerodynamicBeam::AerodynamicBeam");
1509
1510 ASSERT(pBeam != 0);
1511 ASSERT(pBeam->GetElemType() == Elem::BEAM);
1512
1513 for (int iNode = 0; iNode < 3; iNode++) {
1514 pNode[iNode] = pBeam->pGetNode(iNode + 1);
1515
1516 ASSERT(pNode[iNode] != 0);
1517 ASSERT(pNode[iNode]->GetNodeType() == Node::STRUCTURAL);
1518 }
1519 }
1520
~AerodynamicBeam(void)1521 AerodynamicBeam::~AerodynamicBeam(void)
1522 {
1523 DEBUGCOUTFNAME("AerodynamicBeam::~AerodynamicBeam");
1524 }
1525
1526 /* Scrive il contributo dell'elemento al file di restart */
1527 std::ostream&
Restart(std::ostream & out) const1528 AerodynamicBeam::Restart(std::ostream& out) const
1529 {
1530 DEBUGCOUTFNAME("AerodynamicBeam::Restart");
1531 out << " aerodynamic beam: " << GetLabel()
1532 << ", " << pBeam->GetLabel();
1533 if (pIndVel != 0) {
1534 out << ", rotor, " << pIndVel->GetLabel();
1535 }
1536 out << ", reference, node, ", f1.Write(out, ", ")
1537 << ", reference, node, 1, ", (Ra1.GetVec(1)).Write(out, ", ")
1538 << ", 2, ", (Ra1.GetVec(2)).Write(out, ", ")
1539 << ", reference, node, ", f2.Write(out, ", ")
1540 << ", reference, node, 1, ", (Ra2.GetVec(1)).Write(out, ", ")
1541 << ", 2, ", (Ra2.GetVec(2)).Write(out, ", ")
1542 << ", reference, node, ", f3.Write(out, ", ")
1543 << ", reference, node, 1, ", (Ra3.GetVec(1)).Write(out, ", ")
1544 << ", 2, ", (Ra3.GetVec(2)).Write(out, ", ")
1545 << ", ";
1546 Chord.pGetShape()->Restart(out) << ", ";
1547 ForcePoint.pGetShape()->Restart(out) << ", ";
1548 VelocityPoint.pGetShape()->Restart(out) << ", ";
1549 Twist.pGetShape()->Restart(out) << ", "
1550 << GDI.iGetNum() << ", control, ";
1551 pGetDriveCaller()->Restart(out) << ", ";
1552 aerodata->Restart(out);
1553 return out << ";" << std::endl;
1554 }
1555
1556 static const doublereal d13 = 1./sqrt(3.);
1557 static const doublereal pdsi3[] = { -1., -d13, d13 };
1558 static const doublereal pdsf3[] = { -d13, d13, 1. };
1559
1560 /* Jacobian assembly */
1561 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)1562 AerodynamicBeam::AssJac(VariableSubMatrixHandler& WorkMat,
1563 doublereal dCoef,
1564 const VectorHandler& XCurr,
1565 const VectorHandler& XPrimeCurr)
1566 {
1567 DEBUGCOUT("Entering AerodynamicBeam::AssJac()" << std::endl);
1568
1569 if (!bJacobian) {
1570 WorkMat.SetNullMatrix();
1571 return WorkMat;
1572 }
1573
1574 FullSubMatrixHandler& WM = WorkMat.SetFull();
1575
1576 /* Ridimensiona la sottomatrice in base alle esigenze */
1577 integer iNumRows = 0;
1578 integer iNumCols = 0;
1579 WorkSpaceDim(&iNumRows, &iNumCols);
1580 WM.ResizeReset(iNumRows, iNumCols);
1581
1582 integer iNode1FirstIndex = pNode[NODE1]->iGetFirstMomentumIndex();
1583 integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
1584 integer iNode2FirstIndex = pNode[NODE2]->iGetFirstMomentumIndex();
1585 integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
1586 integer iNode3FirstIndex = pNode[NODE3]->iGetFirstMomentumIndex();
1587 integer iNode3FirstPosIndex = pNode[NODE3]->iGetFirstPositionIndex();
1588
1589 for (int iCnt = 1; iCnt <= 6; iCnt++) {
1590 WM.PutRowIndex(iCnt, iNode1FirstIndex + iCnt);
1591 WM.PutColIndex(iCnt, iNode1FirstPosIndex + iCnt);
1592 WM.PutRowIndex(6 + iCnt, iNode2FirstIndex + iCnt);
1593 WM.PutColIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
1594 WM.PutRowIndex(12 + iCnt, iNode3FirstIndex + iCnt);
1595 WM.PutColIndex(12 + iCnt, iNode3FirstPosIndex + iCnt);
1596 }
1597
1598 doublereal dW[6];
1599
1600 /* array di vettori per via del ciclo sui nodi ... */
1601 Vec3 Xn[3];
1602
1603 /* Dati dei nodi */
1604 Xn[NODE1] = pNode[NODE1]->GetXCurr();
1605 const Mat3x3& Rn1(pNode[NODE1]->GetRCurr());
1606 const Vec3& Vn1(pNode[NODE1]->GetVCurr());
1607 const Vec3& Wn1(pNode[NODE1]->GetWCurr());
1608
1609 Xn[NODE2] = pNode[NODE2]->GetXCurr();
1610 const Mat3x3& Rn2(pNode[NODE2]->GetRCurr());
1611 const Vec3& Vn2(pNode[NODE2]->GetVCurr());
1612 const Vec3& Wn2(pNode[NODE2]->GetWCurr());
1613
1614 Xn[NODE3] = pNode[NODE3]->GetXCurr();
1615 const Mat3x3& Rn3(pNode[NODE3]->GetRCurr());
1616 const Vec3& Vn3(pNode[NODE3]->GetVCurr());
1617 const Vec3& Wn3(pNode[NODE3]->GetWCurr());
1618
1619 Vec3 f1Tmp(Rn1*f1);
1620 Vec3 f2Tmp(Rn2*f2);
1621 Vec3 f3Tmp(Rn3*f3);
1622
1623 Vec3 X1Tmp(Xn[NODE1] + f1Tmp);
1624 Vec3 X2Tmp(Xn[NODE2] + f2Tmp);
1625 Vec3 X3Tmp(Xn[NODE3] + f3Tmp);
1626
1627 Vec3 V1Tmp(Vn1 + Wn1.Cross(f1Tmp));
1628 Vec3 V2Tmp(Vn2 + Wn2.Cross(f2Tmp));
1629 Vec3 V3Tmp(Vn3 + Wn3.Cross(f3Tmp));
1630
1631 Vec3 Omega1Crossf1(Wn1.Cross(f1Tmp));
1632 Vec3 Omega2Crossf2(Wn2.Cross(f2Tmp));
1633 Vec3 Omega3Crossf3(Wn3.Cross(f3Tmp));
1634
1635 /*
1636 * Matrice di trasformazione dal sistema globale a quello aerodinamico
1637 */
1638 Mat3x3 RR1(Rn1*Ra1);
1639 Mat3x3 RR2(Rn2*Ra2);
1640 Mat3x3 RR3(Rn3*Ra3);
1641
1642 /*
1643 * Parametri di rotazione dai nodi 1 e 3 al nodo 2 (nell'ipotesi
1644 * che tale trasformazione non dia luogo ad una singolarita')
1645 */
1646
1647 Vec3 g1(ER_Rot::Param, RR2.MulTM(RR1));
1648 Vec3 g3(ER_Rot::Param, RR2.MulTM(RR3));
1649
1650 Mat3x3 GammaInv1(ER_Rot::MatGm1, g1);
1651 Mat3x3 GammaInv3(ER_Rot::MatGm1, g3);
1652
1653 /*
1654 * Se l'elemento e' collegato ad un rotore,
1655 * si fa dare la velocita' di rotazione
1656 */
1657 doublereal dOmega = 0.;
1658 if (pIndVel != 0) {
1659 Rotor *pRotor = dynamic_cast<Rotor *>(pIndVel);
1660 if (pRotor != 0) {
1661 dOmega = pRotor->dGetOmega();
1662 }
1663 }
1664
1665 /*
1666 * Dati "permanenti" (uso solo la posizione del nodo 2 perche'
1667 * non dovrebbero cambiare "molto")
1668 */
1669 doublereal rho, c, p, T;
1670 GetAirProps(Xn[NODE2], rho, c, p, T); /* p, T no used yet */
1671 aerodata->SetAirData(rho, c);
1672
1673 int iPnt = 0;
1674
1675 ResetIterator();
1676
1677 integer iNumDof = aerodata->iGetNumDof();
1678 integer iFirstEq = -1;
1679 integer iFirstSubEq = -1;
1680 if (iNumDof > 0) {
1681 iFirstEq = iGetFirstIndex();
1682 iFirstSubEq = 18;
1683
1684 integer iOffset = iFirstEq - 18;
1685
1686 for (int iCnt = 18 + 1; iCnt <= iNumRows; iCnt++) {
1687 WM.PutRowIndex(iCnt, iOffset + iCnt);
1688 WM.PutColIndex(iCnt, iOffset + iCnt);
1689 }
1690 }
1691
1692 for (int iNode = 0; iNode < LASTNODE; iNode++) {
1693 doublereal dsi = pdsi3[iNode];
1694 doublereal dsf = pdsf3[iNode];
1695
1696 doublereal dsm = (dsf + dsi)/2.;
1697 doublereal dsdCsi = (dsf - dsi)/2.;
1698
1699 Mat3x3 WM_F[6], WM_M[6];
1700 WM_F[DELTAx1].Reset();
1701 WM_F[DELTAg1].Reset();
1702 WM_F[DELTAx2].Reset();
1703 WM_F[DELTAg2].Reset();
1704 WM_F[DELTAx3].Reset();
1705 WM_F[DELTAg3].Reset();
1706 WM_M[DELTAx1].Reset();
1707 WM_M[DELTAg1].Reset();
1708 WM_M[DELTAx2].Reset();
1709 WM_M[DELTAg2].Reset();
1710 WM_M[DELTAx3].Reset();
1711 WM_M[DELTAg3].Reset();
1712
1713 //unsigned int iDelta_x1, iDelta_g1, iDelta_x2, iDelta_g2, iDelta_x3, iDelta_g3;
1714 //iDelta_x1 = 0;, iDelta_g1, iDelta_x2, iDelta_g2, iDelta_x3, iDelta_g3;
1715
1716 /* Ciclo sui punti di Gauss */
1717 PntWght PW = GDI.GetFirst();
1718 do {
1719 doublereal dCsi = PW.dGetPnt();
1720 doublereal ds = dsm + dsdCsi*dCsi;
1721 doublereal dXds = DxDcsi3N(ds,
1722 Xn[NODE1], Xn[NODE2], Xn[NODE3]);
1723
1724 doublereal dN1 = ShapeFunc3N(ds, 1);
1725 doublereal dN2 = ShapeFunc3N(ds, 2);
1726 doublereal dN3 = ShapeFunc3N(ds, 3);
1727
1728 Vec3 Xr(X1Tmp*dN1 + X2Tmp*dN2 + X3Tmp*dN3);
1729 Vec3 Vr(V1Tmp*dN1 + V2Tmp*dN2 + V3Tmp*dN3);
1730 Vec3 Wr(Wn1*dN1 + Wn2*dN2 + Wn3*dN3);
1731 Vec3 gr(g1*dN1 + g3*dN3);
1732 Mat3x3 Gamma(ER_Rot::MatG, gr);
1733
1734 /* Contributo di velocita' del vento */
1735 Vec3 VTmp(Zero3);
1736 if (fGetAirVelocity(VTmp, Xr)) {
1737 Vr -= VTmp;
1738 }
1739
1740 /*
1741 * Se l'elemento e' collegato ad un rotore,
1742 * aggiunge alla velocita' la velocita' indotta
1743 */
1744 if (pIndVel != 0) {
1745 Vr += pIndVel->GetInducedVelocity(GetElemType(),
1746 GetLabel(), iPnt, Xr);
1747 }
1748
1749 /* Copia i dati nel vettore di lavoro dVAM */
1750 doublereal dTw = Twist.dGet(ds);
1751 /* Contributo dell'eventuale sup. mobile */
1752 dTw += dGet();
1753
1754 aerodata->SetSectionData(dCsi,
1755 Chord.dGet(ds),
1756 ForcePoint.dGet(ds),
1757 VelocityPoint.dGet(ds),
1758 dTw,
1759 dOmega);
1760
1761 /*
1762 * Lo svergolamento non viene piu' trattato in aerod2_;
1763 * quindi lo uso per correggere la matrice di rotazione
1764 * dal sistema aerodinamico a quello globale
1765 */
1766 Mat3x3 RRloc(RR2*Mat3x3(ER_Rot::MatR, gr));
1767 if (dTw != 0.) {
1768 doublereal dCosT = cos(dTw);
1769 doublereal dSinT = sin(dTw);
1770 /* Assumo lo svergolamento positivo a cabrare */
1771 Mat3x3 RTw(dCosT, dSinT, 0.,
1772 -dSinT, dCosT, 0.,
1773 0., 0., 1.);
1774 /*
1775 * Allo stesso tempo interpola le g
1776 * e aggiunge lo svergolamento
1777 */
1778 RRloc = RRloc*RTw;
1779 }
1780
1781 /*
1782 * Ruota velocita' e velocita' angolare nel sistema
1783 * aerodinamico e li copia nel vettore di lavoro dW
1784 */
1785 VTmp = RRloc.MulTV(Vr);
1786 VTmp.PutTo(&dW[0]);
1787
1788 Vec3 WTmp = RRloc.MulTV(Wr);
1789 WTmp.PutTo(&dW[3]);
1790 /* Funzione di calcolo delle forze aerodinamiche */
1791 doublereal Fa0[6];
1792 Mat6x6 JFa;
1793
1794 doublereal cc = dXds*dsdCsi*PW.dGetWght();
1795
1796 Vec3 d(Xr - Xn[iNode]);
1797
1798 Mat3x3 Theta1(RR2*Gamma*GammaInv1.MulMT(RR2*dN1));
1799 Mat3x3 Theta3(RR2*Gamma*GammaInv3.MulMT(RR2*dN3));
1800 Mat3x3 Theta2(Eye3 - Theta1 - Theta3);
1801
1802 Vec3 Vrc(Vr*dCoef);
1803 Mat3x3 Bv1(Vrc.Cross(Theta1) - Mat3x3(MatCross, Omega1Crossf1*(dN1*dCoef)));
1804 Mat3x3 Bv2(Vrc.Cross(Theta2) - Mat3x3(MatCross, Omega2Crossf2*(dN2*dCoef)));
1805 Mat3x3 Bv3(Vrc.Cross(Theta3) - Mat3x3(MatCross, Omega3Crossf3*(dN3*dCoef)));
1806
1807 Vec3 Wrc(Wr*dCoef);
1808 Mat3x3 Bw1(Wrc.Cross(Theta1) - Mat3x3(MatCross, Wn1*(dN1*dCoef)));
1809 Mat3x3 Bw2(Wrc.Cross(Theta2) - Mat3x3(MatCross, Wn2*(dN2*dCoef)));
1810 Mat3x3 Bw3(Wrc.Cross(Theta3) - Mat3x3(MatCross, Wn3*(dN3*dCoef)));
1811
1812 if (iNumDof) {
1813 // prepare (v/dot{x} + dCoef*v/x) and so
1814 Mat3x3 RRlocT(RRloc.Transpose());
1815
1816 vx.PutMat3x3(1, RRlocT*dN1);
1817 vx.PutMat3x3(4, RRloc.MulTM(Bv1 - Mat3x3(MatCross, f1Tmp*dN1)));
1818
1819 vx.PutMat3x3(6 + 1, RRlocT*dN2);
1820 vx.PutMat3x3(6 + 4, RRloc.MulTM(Bv2 - Mat3x3(MatCross, f2Tmp*dN2)));
1821
1822 vx.PutMat3x3(12 + 1, RRlocT*dN3);
1823 vx.PutMat3x3(12 + 4, RRloc.MulTM(Bv3 - Mat3x3(MatCross, f3Tmp*dN3)));
1824
1825 wx.PutMat3x3(4, RRlocT + Bw1);
1826 wx.PutMat3x3(6 + 4, RRlocT + Bw2);
1827 wx.PutMat3x3(12 + 4, RRlocT + Bw3);
1828
1829 // equations from iFirstEq on are dealt with by aerodata
1830 aerodata->AssJac(WM, dCoef, XCurr, XPrimeCurr,
1831 iFirstEq, iFirstSubEq,
1832 vx, wx, fq, cq, iPnt, dW, Fa0, JFa, OUTA[iPnt]);
1833
1834 // deal with (f/dot{q} + dCoef*f/q) and so
1835 integer iOffset = 18 + iPnt*iNumDof;
1836 for (integer iCol = 1; iCol <= iNumDof; iCol++) {
1837 Vec3 fqTmp((RRloc*fq.GetVec(iCol))*cc);
1838 Vec3 cqTmp(d.Cross(fqTmp) + (RRloc*cq.GetVec(iCol))*cc);
1839
1840 WM.Sub(6*iNode + 1, iOffset + iCol, fqTmp);
1841 WM.Sub(6*iNode + 4, iOffset + iCol, cqTmp);
1842 }
1843
1844 // first equation
1845 iFirstEq += iNumDof;
1846 iFirstSubEq += iNumDof;
1847
1848 } else {
1849 aerodata->GetForcesJac(iPnt, dW, Fa0, JFa, OUTA[iPnt]);
1850 }
1851
1852 // rotate force, couple and Jacobian matrix in absolute frame
1853 Mat6x6 JFaR = MultRMRt(JFa, RRloc, cc);
1854
1855 // force and moment about the node
1856 Vec3 fTmp(RRloc*(Vec3(&Fa0[0])*dCoef));
1857 Vec3 cTmp(RRloc*(Vec3(&Fa0[3])*dCoef) + d.Cross(fTmp));
1858
1859 Mat3x3 WM_F2[6];
1860
1861 // f <-> x
1862 WM_F2[DELTAx1] = JFaR.GetMat11()*dN1;
1863
1864 WM_F2[DELTAx2] = JFaR.GetMat11()*dN2;
1865
1866 WM_F2[DELTAx3] = JFaR.GetMat11()*dN3;
1867
1868 doublereal delta;
1869
1870 // c <-> x
1871 delta = (iNode == NODE1) ? 1. : 0.;
1872 WM_M[DELTAx1] += JFaR.GetMat21()*dN1 - Mat3x3(MatCross, fTmp*(dN1 - delta));
1873
1874 delta = (iNode == NODE2) ? 1. : 0.;
1875 WM_M[DELTAx2] += JFaR.GetMat21()*dN2 - Mat3x3(MatCross, fTmp*(dN2 - delta));
1876
1877 delta = (iNode == NODE3) ? 1. : 0.;
1878 WM_M[DELTAx3] += JFaR.GetMat21()*dN3 - Mat3x3(MatCross, fTmp*(dN3 - delta));
1879
1880 // f <-> g
1881 WM_F2[DELTAg1] = (JFaR.GetMat12() - JFaR.GetMat11()*Mat3x3(MatCross, f1Tmp))*dN1;
1882 WM_F2[DELTAg1] += JFaR.GetMat11()*Bv1 + JFaR.GetMat12()*Bw1;
1883 WM_F2[DELTAg1] -= fTmp.Cross(Theta1);
1884
1885 WM_F2[DELTAg2] = (JFaR.GetMat12() - JFaR.GetMat11()*Mat3x3(MatCross, f2Tmp))*dN2;
1886 WM_F2[DELTAg2] += JFaR.GetMat11()*Bv2 + JFaR.GetMat12()*Bw2;
1887 WM_F2[DELTAg2] -= fTmp.Cross(Theta2);
1888
1889 WM_F2[DELTAg3] = (JFaR.GetMat12() - JFaR.GetMat11()*Mat3x3(MatCross, f3Tmp))*dN3;
1890 WM_F2[DELTAg3] += JFaR.GetMat11()*Bv3 + JFaR.GetMat12()*Bw3;
1891 WM_F2[DELTAg3] -= fTmp.Cross(Theta3);
1892
1893 // c <-> g
1894 WM_M[DELTAg1] += (JFaR.GetMat22() - JFaR.GetMat21()*Mat3x3(MatCross, f1Tmp))*dN1;
1895 WM_M[DELTAg1] += JFaR.GetMat21()*Bv1 + JFaR.GetMat22()*Bw1;
1896 WM_M[DELTAg1] -= cTmp.Cross(Theta1);
1897 WM_M[DELTAg1] += Mat3x3(MatCrossCross, fTmp, f1Tmp*dN1);
1898
1899 WM_M[DELTAg2] += (JFaR.GetMat22() - JFaR.GetMat21()*Mat3x3(MatCross, f2Tmp))*dN2;
1900 WM_M[DELTAg2] += JFaR.GetMat21()*Bv2 + JFaR.GetMat22()*Bw2;
1901 WM_M[DELTAg2] -= cTmp.Cross(Theta2);
1902 WM_M[DELTAg2] += Mat3x3(MatCrossCross, fTmp, f2Tmp*dN2);
1903
1904 WM_M[DELTAg3] += (JFaR.GetMat22() - JFaR.GetMat21()*Mat3x3(MatCross, f3Tmp))*dN3;
1905 WM_M[DELTAg3] += JFaR.GetMat21()*Bv3 + JFaR.GetMat22()*Bw3;
1906 WM_M[DELTAg3] -= cTmp.Cross(Theta3);
1907 WM_M[DELTAg3] += Mat3x3(MatCrossCross, fTmp, f3Tmp*dN3);
1908
1909 for (int iCnt = 0; iCnt < 2*LASTNODE; iCnt++) {
1910 WM_F[iCnt] += WM_F2[iCnt];
1911 WM_M[iCnt] += d.Cross(WM_F2[iCnt]);
1912 }
1913
1914 iPnt++;
1915
1916 } while (GDI.fGetNext(PW));
1917
1918 for (int iCnt = 0; iCnt < 2*LASTNODE; iCnt++) {
1919 WM.Sub(6*iNode + 1, 3*iCnt + 1, WM_F[iCnt]);
1920 WM.Sub(6*iNode + 4, 3*iCnt + 1, WM_M[iCnt]);
1921 }
1922 }
1923
1924 return WorkMat;
1925 }
1926
1927 /* assemblaggio residuo */
1928 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)1929 AerodynamicBeam::AssRes(SubVectorHandler& WorkVec,
1930 doublereal dCoef,
1931 const VectorHandler& XCurr,
1932 const VectorHandler& XPrimeCurr)
1933 {
1934 DEBUGCOUTFNAME("AerodynamicBeam::AssRes");
1935
1936 integer iNumRows;
1937 integer iNumCols;
1938 WorkSpaceDim(&iNumRows, &iNumCols);
1939 WorkVec.ResizeReset(iNumRows);
1940
1941 integer iNode1FirstIndex = pNode[NODE1]->iGetFirstMomentumIndex();
1942 integer iNode2FirstIndex = pNode[NODE2]->iGetFirstMomentumIndex();
1943 integer iNode3FirstIndex = pNode[NODE3]->iGetFirstMomentumIndex();
1944 for (int iCnt = 1; iCnt <= 6; iCnt++) {
1945 WorkVec.PutRowIndex(iCnt, iNode1FirstIndex + iCnt);
1946 WorkVec.PutRowIndex(6 + iCnt, iNode2FirstIndex + iCnt);
1947 WorkVec.PutRowIndex(12 + iCnt, iNode3FirstIndex + iCnt);
1948 }
1949
1950 AssVec(WorkVec, dCoef, XCurr, XPrimeCurr);
1951
1952 return WorkVec;
1953 }
1954
1955 /* assemblaggio iniziale residuo */
1956 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)1957 AerodynamicBeam::InitialAssRes(SubVectorHandler& WorkVec,
1958 const VectorHandler& XCurr)
1959 {
1960 DEBUGCOUTFNAME("AerodynamicBeam::InitialAssRes");
1961 WorkVec.ResizeReset(18);
1962
1963 integer iNode1FirstIndex = pNode[NODE1]->iGetFirstPositionIndex();
1964 integer iNode2FirstIndex = pNode[NODE2]->iGetFirstPositionIndex();
1965 integer iNode3FirstIndex = pNode[NODE3]->iGetFirstPositionIndex();
1966 for (int iCnt = 1; iCnt <= 6; iCnt++) {
1967 WorkVec.PutRowIndex(iCnt, iNode1FirstIndex + iCnt);
1968 WorkVec.PutRowIndex(6 + iCnt, iNode2FirstIndex + iCnt);
1969 WorkVec.PutRowIndex(12 + iCnt, iNode3FirstIndex + iCnt);
1970 }
1971
1972 AssVec(WorkVec, 1., XCurr, XCurr);
1973
1974 return WorkVec;
1975 }
1976
1977
1978 /* assemblaggio residuo */
1979 void
AssVec(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)1980 AerodynamicBeam::AssVec(SubVectorHandler& WorkVec,
1981 doublereal dCoef,
1982 const VectorHandler& XCurr,
1983 const VectorHandler& XPrimeCurr)
1984 {
1985 DEBUGCOUTFNAME("AerodynamicBeam::AssVec");
1986
1987 /* array di vettori per via del ciclo sui nodi ... */
1988 Vec3 Xn[3];
1989
1990 /* Dati dei nodi */
1991 Xn[NODE1] = pNode[NODE1]->GetXCurr();
1992 const Mat3x3& Rn1(pNode[NODE1]->GetRCurr());
1993 const Vec3& Vn1(pNode[NODE1]->GetVCurr());
1994 const Vec3& Wn1(pNode[NODE1]->GetWCurr());
1995
1996 Xn[NODE2] = pNode[NODE2]->GetXCurr();
1997 const Mat3x3& Rn2(pNode[NODE2]->GetRCurr());
1998 const Vec3& Vn2(pNode[NODE2]->GetVCurr());
1999 const Vec3& Wn2(pNode[NODE2]->GetWCurr());
2000
2001 Xn[NODE3] = pNode[NODE3]->GetXCurr();
2002 const Mat3x3& Rn3(pNode[NODE3]->GetRCurr());
2003 const Vec3& Vn3(pNode[NODE3]->GetVCurr());
2004 const Vec3& Wn3(pNode[NODE3]->GetWCurr());
2005
2006 Vec3 f1Tmp(Rn1*f1);
2007 Vec3 f2Tmp(Rn2*f2);
2008 Vec3 f3Tmp(Rn3*f3);
2009
2010 Vec3 X1Tmp(Xn[NODE1] + f1Tmp);
2011 Vec3 X2Tmp(Xn[NODE2] + f2Tmp);
2012 Vec3 X3Tmp(Xn[NODE3] + f3Tmp);
2013
2014 Vec3 V1Tmp(Vn1 + Wn1.Cross(f1Tmp));
2015 Vec3 V2Tmp(Vn2 + Wn2.Cross(f2Tmp));
2016 Vec3 V3Tmp(Vn3 + Wn3.Cross(f3Tmp));
2017
2018 /*
2019 * Matrice di trasformazione dal sistema globale a quello aerodinamico
2020 */
2021 Mat3x3 RR1(Rn1*Ra1);
2022 Mat3x3 RR2(Rn2*Ra2);
2023 Mat3x3 RR3(Rn3*Ra3);
2024
2025 /*
2026 * Parametri di rotazione dai nodi 1 e 3 al nodo 2 (nell'ipotesi
2027 * che tale trasformazione non dia luogo ad una singolarita')
2028 */
2029 Vec3 g1(ER_Rot::Param, RR2.MulTM(RR1));
2030 Vec3 g3(ER_Rot::Param, RR2.MulTM(RR3));
2031
2032 /*
2033 * Se l'elemento e' collegato ad un rotore,
2034 * si fa dare la velocita' di rotazione
2035 */
2036 doublereal dOmega = 0.;
2037 if (pIndVel != 0) {
2038 Rotor *pRotor = dynamic_cast<Rotor *>(pIndVel);
2039 if (pRotor != 0) {
2040 dOmega = pRotor->dGetOmega();
2041 }
2042 }
2043
2044 /*
2045 * Dati "permanenti" (uso solo la posizione del nodo 2 perche'
2046 * non dovrebbero cambiare "molto")
2047 */
2048 doublereal rho, c, p, T;
2049 GetAirProps(Xn[NODE2], rho, c, p, T); /* p, T no used yet */
2050 aerodata->SetAirData(rho, c);
2051
2052 int iPnt = 0;
2053
2054 ResetIterator();
2055
2056 integer iNumDof = aerodata->iGetNumDof();
2057 integer iFirstEq = -1;
2058 integer iFirstSubEq = -1;
2059 if (iNumDof > 0) {
2060 iFirstEq = iGetFirstIndex();
2061 iFirstSubEq = 18;
2062
2063 integer iOffset = iFirstEq - 18;
2064 integer iNumRows = WorkVec.iGetSize();
2065 for (int iCnt = 18 + 1; iCnt <= iNumRows; iCnt++) {
2066 WorkVec.PutRowIndex(iCnt, iOffset + iCnt);
2067 }
2068 }
2069
2070 for (int iNode = 0; iNode < LASTNODE; iNode++) {
2071
2072 /* Resetta le forze */
2073 F[iNode].Reset();
2074 M[iNode].Reset();
2075
2076 doublereal dsi = pdsi3[iNode];
2077 doublereal dsf = pdsf3[iNode];
2078
2079 doublereal dsm = (dsf + dsi)/2.;
2080 doublereal dsdCsi = (dsf - dsi)/2.;
2081
2082 /* Ciclo sui punti di Gauss */
2083 PntWght PW = GDI.GetFirst();
2084 do {
2085 doublereal dCsi = PW.dGetPnt();
2086 doublereal ds = dsm + dsdCsi*dCsi;
2087 doublereal dXds = DxDcsi3N(ds,
2088 Xn[NODE1], Xn[NODE2], Xn[NODE3]);
2089
2090 doublereal dN1 = ShapeFunc3N(ds, 1);
2091 doublereal dN2 = ShapeFunc3N(ds, 2);
2092 doublereal dN3 = ShapeFunc3N(ds, 3);
2093
2094 Vec3 Xr(X1Tmp*dN1 + X2Tmp*dN2 + X3Tmp*dN3);
2095 Vec3 Vr(V1Tmp*dN1 + V2Tmp*dN2 + V3Tmp*dN3);
2096 Vec3 Wr(Wn1*dN1 + Wn2*dN2 + Wn3*dN3);
2097
2098 /* Contributo di velocita' del vento */
2099 Vec3 VTmp(Zero3);
2100 if (fGetAirVelocity(VTmp, Xr)) {
2101 Vr -= VTmp;
2102 }
2103
2104 /*
2105 * Se l'elemento e' collegato ad un rotore,
2106 * aggiunge alla velocita' la velocita' indotta
2107 */
2108 if (pIndVel != 0) {
2109 Vr += pIndVel->GetInducedVelocity(GetElemType(),
2110 GetLabel(), iPnt, Xr);
2111 }
2112
2113 /* Copia i dati nel vettore di lavoro dVAM */
2114 doublereal dTw = Twist.dGet(ds);
2115 /* Contributo dell'eventuale sup. mobile */
2116 dTw += dGet();
2117
2118 aerodata->SetSectionData(dCsi,
2119 Chord.dGet(ds),
2120 ForcePoint.dGet(ds),
2121 VelocityPoint.dGet(ds),
2122 dTw,
2123 dOmega);
2124
2125 /*
2126 * Lo svergolamento non viene piu' trattato in aerod2_;
2127 * quindi lo uso per correggere la matrice di rotazione
2128 * dal sistema aerodinamico a quello globale
2129 */
2130 Mat3x3 RRloc(RR2*Mat3x3(ER_Rot::MatR, g1*dN1 + g3*dN3));
2131 if (dTw != 0.) {
2132 doublereal dCosT = cos(dTw);
2133 doublereal dSinT = sin(dTw);
2134 /* Assumo lo svergolamento positivo a cabrare */
2135 Mat3x3 RTw(dCosT, dSinT, 0.,
2136 -dSinT, dCosT, 0.,
2137 0., 0., 1.);
2138 /*
2139 * Allo stesso tempo interpola le g
2140 * e aggiunge lo svergolamento
2141 */
2142 RRloc = RRloc*RTw;
2143 }
2144
2145 /*
2146 * Ruota velocita' e velocita' angolare nel sistema
2147 * aerodinamico e li copia nel vettore di lavoro dW
2148 */
2149 doublereal dW[6];
2150 doublereal dTng[6];
2151
2152 VTmp = RRloc.MulTV(Vr);
2153 VTmp.PutTo(&dW[0]);
2154
2155 Vec3 WTmp = RRloc.MulTV(Wr);
2156 WTmp.PutTo(&dW[3]);
2157
2158 /* Funzione di calcolo delle forze aerodinamiche */
2159 if (iNumDof) {
2160 aerodata->AssRes(WorkVec, dCoef, XCurr, XPrimeCurr,
2161 iFirstEq, iFirstSubEq, iPnt, dW, dTng, OUTA[iPnt]);
2162
2163 // first equation
2164 iFirstEq += iNumDof;
2165 iFirstSubEq += iNumDof;
2166
2167 } else {
2168 aerodata->GetForces(iPnt, dW, dTng, OUTA[iPnt]);
2169 }
2170
2171 /* Dimensionalizza le forze */
2172 doublereal dWght = dXds*dsdCsi*PW.dGetWght();
2173 dTng[1] *= TipLoss.dGet(dCsi);
2174 Vec3 FTmp(RRloc*(Vec3(&dTng[0])));
2175 Vec3 MTmp(RRloc*(Vec3(&dTng[3])));
2176
2177 // Se e' definito il rotore, aggiungere il contributo alla trazione
2178 AddSectionalForce_int(iPnt, FTmp, MTmp, dWght, Xr, RRloc, Vr, Wr);
2179
2180 FTmp *= dWght;
2181 MTmp *= dWght;
2182 F[iNode] += FTmp;
2183 M[iNode] += MTmp;
2184 M[iNode] += (Xr - Xn[iNode]).Cross(FTmp);
2185
2186 // specific for Gauss points force output
2187 if (bToBeOutput()) {
2188 SetData(VTmp, dTng, Xr, RRloc, Vr, Wr, FTmp, MTmp);
2189 }
2190
2191 iPnt++;
2192
2193 } while (GDI.fGetNext(PW));
2194
2195 // Se e' definito il rotore, aggiungere il contributo alla trazione
2196 AddForce_int(pNode[iNode], F[iNode], M[iNode], Xn[iNode]);
2197
2198 /* Somma il termine al residuo */
2199 WorkVec.Add(6*iNode + 1, F[iNode]);
2200 WorkVec.Add(6*iNode + 4, M[iNode]);
2201 }
2202 }
2203
2204 /*
2205 * output; si assume che ogni tipo di elemento sappia, attraverso
2206 * l'OutputHandler, dove scrivere il proprio output
2207 */
2208 void
Output(OutputHandler & OH) const2209 AerodynamicBeam::Output(OutputHandler& OH) const
2210 {
2211 DEBUGCOUTFNAME("AerodynamicBeam::Output");
2212
2213 if (bToBeOutput()) {
2214 Aerodynamic2DElem<3>::Output_int(OH);
2215
2216 if (OH.UseText(OutputHandler::AERODYNAMIC)) {
2217 std::ostream& out = OH.Aerodynamic() << std::setw(8) << GetLabel();
2218
2219 switch (GetOutput()) {
2220 case AEROD_OUT_NODE:
2221 out << " " << std::setw(8) << pBeam->GetLabel()
2222 << " ", F[NODE1].Write(out, " ") << " ", M[NODE1].Write(out, " ")
2223 << " ", F[NODE2].Write(out, " ") << " ", M[NODE2].Write(out, " ")
2224 << " ", F[NODE3].Write(out, " ") << " ", M[NODE3].Write(out, " ");
2225 break;
2226
2227 case AEROD_OUT_PGAUSS:
2228 ASSERT(!OutputData.empty());
2229
2230 for (std::vector<Aero_output>::const_iterator i = OutputData.begin();
2231 i != OutputData.end(); ++i)
2232 {
2233 out << " " << i->alpha
2234 << " " << i->f;
2235 }
2236 break;
2237
2238 case AEROD_OUT_STD:
2239 for (int i = 0; i < 3*GDI.iGetNum(); i++) {
2240 out
2241 << " " << OUTA[i].alpha
2242 << " " << OUTA[i].gamma
2243 << " " << OUTA[i].mach
2244 << " " << OUTA[i].cl
2245 << " " << OUTA[i].cd
2246 << " " << OUTA[i].cm
2247 << " " << OUTA[i].alf1
2248 << " " << OUTA[i].alf2;
2249 }
2250 break;
2251
2252 default:
2253 ASSERT(0);
2254 break;
2255 }
2256
2257 out << std::endl;
2258 }
2259 }
2260 }
2261
2262 /* AerodynamicBeam - end */
2263
2264
2265 /* Legge un elemento aerodinamico di trave */
2266
2267 Elem *
ReadAerodynamicBeam(DataManager * pDM,MBDynParser & HP,const DofOwner * pDO,unsigned int uLabel)2268 ReadAerodynamicBeam(DataManager* pDM,
2269 MBDynParser& HP,
2270 const DofOwner *pDO,
2271 unsigned int uLabel)
2272 {
2273 DEBUGCOUTFNAME("ReadAerodynamicBeam");
2274
2275 /* Trave */
2276 unsigned int uBeam = (unsigned int)HP.GetInt();
2277
2278 DEBUGLCOUT(MYDEBUG_INPUT, "Linked to beam: " << uBeam << std::endl);
2279
2280 /* verifica di esistenza della trave */
2281 Elem* p = pDM->pFindElem(Elem::BEAM, uBeam);
2282 if (p == 0) {
2283 silent_cerr("Beam3(" << uBeam << ") not defined "
2284 "at line " << HP.GetLineData()
2285 << std::endl);
2286 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2287 }
2288 Beam *pBeam = dynamic_cast<Beam *>(p);
2289 if (pBeam == 0) {
2290 silent_cerr("Beam(" << uBeam << ") is not a Beam3 "
2291 "at line " << HP.GetLineData()
2292 << std::endl);
2293 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2294 }
2295 ASSERT(pBeam != 0);
2296
2297 /* Eventuale rotore */
2298 InducedVelocity* pIndVel = 0;
2299 bool bPassive(false);
2300 (void)ReadInducedVelocity(pDM, HP, uLabel, "AerodynamicBeam3",
2301 pIndVel, bPassive);
2302
2303 /* Nodo 1: */
2304
2305 /* Offset del corpo aerodinamico rispetto al nodo */
2306 const StructNode* pNode1 = pBeam->pGetNode(1);
2307
2308 ReferenceFrame RF(pNode1);
2309 Vec3 f1(HP.GetPosRel(RF));
2310 DEBUGLCOUT(MYDEBUG_INPUT, "Node 1 offset: " << f1 << std::endl);
2311
2312 Mat3x3 Ra1(HP.GetRotRel(RF));
2313 DEBUGLCOUT(MYDEBUG_INPUT,
2314 "Node 1 rotation matrix: " << std::endl << Ra1 << std::endl);
2315
2316 /* Nodo 2: */
2317
2318 /* Offset del corpo aerodinamico rispetto al nodo */
2319 const StructNode* pNode2 = pBeam->pGetNode(2);
2320
2321 RF = ReferenceFrame(pNode2);
2322 Vec3 f2(HP.GetPosRel(RF));
2323 DEBUGLCOUT(MYDEBUG_INPUT, "Node 2 offset: " << f2 << std::endl);
2324
2325 Mat3x3 Ra2(HP.GetRotRel(RF));
2326 DEBUGLCOUT(MYDEBUG_INPUT,
2327 "Node 2 rotation matrix: " << std::endl << Ra2 << std::endl);
2328
2329 /* Nodo 3: */
2330
2331 /* Offset del corpo aerodinamico rispetto al nodo */
2332 const StructNode* pNode3 = pBeam->pGetNode(3);
2333
2334 RF = ReferenceFrame(pNode3);
2335 Vec3 f3(HP.GetPosRel(RF));
2336 DEBUGLCOUT(MYDEBUG_INPUT, "Node 3 offset: " << f3 << std::endl);
2337
2338 Mat3x3 Ra3(HP.GetRotRel(RF));
2339 DEBUGLCOUT(MYDEBUG_INPUT,
2340 "Node 3 rotation matrix: " << std::endl << Ra3 << std::endl);
2341
2342 Shape* pChord = 0;
2343 Shape* pForce = 0;
2344 Shape* pVelocity = 0;
2345 Shape* pTwist = 0;
2346 Shape* pTipLoss = 0;
2347
2348 integer iNumber = 0;
2349 DriveCaller* pDC = 0;
2350 AeroData* aerodata = 0;
2351
2352 ReadAeroData(pDM, HP, 3,
2353 &pChord, &pForce, &pVelocity, &pTwist, &pTipLoss,
2354 &iNumber, &pDC, &aerodata);
2355
2356 bool bUseJacobian(false);
2357 if (HP.IsKeyWord("jacobian")) {
2358 bUseJacobian = HP.GetYesNoOrBool(bDefaultUseJacobian);
2359 }
2360
2361 if (aerodata->iGetNumDof() > 0 && !bUseJacobian) {
2362 silent_cerr("AerodynamicBeam3(" << uLabel << "): "
2363 "aerodynamic model needs \"jacobian, yes\" at line " << HP.GetLineData()
2364 << std::endl);
2365 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2366 }
2367
2368 OrientationDescription od = UNKNOWN_ORIENTATION_DESCRIPTION;
2369 unsigned uFlags = AerodynamicOutput::OUTPUT_NONE;
2370 ReadOptionalAerodynamicCustomOutput(pDM, HP, uLabel, uFlags, od);
2371
2372 flag fOut = pDM->fReadOutput(HP, Elem::AERODYNAMIC);
2373 if (HP.IsArg()) {
2374 if (HP.IsKeyWord("std")) {
2375 fOut |= AerodynamicOutput::AEROD_OUT_STD;
2376 } else if (HP.IsKeyWord("gauss")) {
2377 fOut |= AerodynamicOutput::AEROD_OUT_PGAUSS;
2378 } else if (HP.IsKeyWord("node")) {
2379 fOut |= AerodynamicOutput::AEROD_OUT_NODE;
2380 } else {
2381 silent_cerr("AerodynamicBeam3(" << uLabel << "): "
2382 "unknown output mode at line " << HP.GetLineData()
2383 << std::endl);
2384 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2385 }
2386
2387 } else if (fOut) {
2388 fOut |= AerodynamicOutput::AEROD_OUT_STD;
2389 }
2390 fOut |= uFlags;
2391
2392 Elem* pEl = 0;
2393
2394 SAFENEWWITHCONSTRUCTOR(pEl,
2395 AerodynamicBeam,
2396 AerodynamicBeam(uLabel, pDO, pBeam, pIndVel, bPassive,
2397 f1, f2, f3, Ra1, Ra2, Ra3,
2398 pChord, pForce, pVelocity, pTwist, pTipLoss,
2399 iNumber, aerodata, pDC, bUseJacobian, od, fOut));
2400
2401 /* Se non c'e' il punto e virgola finale */
2402 if (HP.IsArg()) {
2403 silent_cerr("semicolon expected at line "
2404 << HP.GetLineData() << std::endl);
2405 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2406 }
2407
2408 std::ostream& out = pDM->GetLogFile();
2409 out << "aero3: " << uLabel;
2410
2411 Vec3 ra1 = Ra1.GetVec(1);
2412 Vec3 ra3 = Ra1.GetVec(3);
2413 doublereal dC = pChord->dGet(-1.);
2414 doublereal dP = pForce->dGet(-1.);
2415 out
2416 << " " << pNode1->GetLabel()
2417 << " ", (f1 + ra1*(dP - dC*3./4.)).Write(out, " ")
2418 << " ", (f1 + ra1*(dP + dC/4.)).Write(out, " ");
2419
2420 ra1 = Ra2.GetVec(1);
2421 ra3 = Ra2.GetVec(3);
2422 dC = pChord->dGet(0.);
2423 dP = pForce->dGet(0.);
2424 out
2425 << " " << pNode2->GetLabel()
2426 << " ", (f2 + ra1*(dP - dC*3./4.)).Write(out, " ")
2427 << " ", (f2 + ra1*(dP + dC/4.)).Write(out, " ");
2428
2429 ra1 = Ra3.GetVec(1);
2430 ra3 = Ra3.GetVec(3);
2431 dC = pChord->dGet(1.);
2432 dP = pForce->dGet(1.);
2433 out
2434 << " " << pNode3->GetLabel()
2435 << " ", (f3 + ra1*(dP - dC*3./4.)).Write(out, " ")
2436 << " ", (f3 + ra1*(dP + dC/4.)).Write(out, " ")
2437 << std::endl;
2438
2439 return pEl;
2440 } /* End of ReadAerodynamicBeam() */
2441
2442
2443 /* AerodynamicBeam2 - begin */
2444
AerodynamicBeam2(unsigned int uLabel,const DofOwner * pDO,const Beam2 * pB,InducedVelocity * pR,bool bPassive,const Vec3 & fTmp1,const Vec3 & fTmp2,const Mat3x3 & Ra1Tmp,const Mat3x3 & Ra2Tmp,const Shape * pC,const Shape * pF,const Shape * pV,const Shape * pT,const Shape * pTL,integer iN,AeroData * a,const DriveCaller * pDC,bool bUseJacobian,OrientationDescription ood,flag fOut)2445 AerodynamicBeam2::AerodynamicBeam2(
2446 unsigned int uLabel,
2447 const DofOwner* pDO,
2448 const Beam2* pB,
2449 InducedVelocity* pR, bool bPassive,
2450 const Vec3& fTmp1,
2451 const Vec3& fTmp2,
2452 const Mat3x3& Ra1Tmp,
2453 const Mat3x3& Ra2Tmp,
2454 const Shape* pC,
2455 const Shape* pF,
2456 const Shape* pV,
2457 const Shape* pT,
2458 const Shape* pTL,
2459 integer iN,
2460 AeroData* a,
2461 const DriveCaller* pDC,
2462 bool bUseJacobian,
2463 OrientationDescription ood,
2464 flag fOut
2465 )
2466 : Elem(uLabel, fOut),
2467 Aerodynamic2DElem<2>(uLabel, pDO, pR, bPassive,
2468 pC, pF, pV, pT, pTL, iN, a, pDC, bUseJacobian, ood, fOut),
2469 pBeam(pB),
2470 f1(fTmp1),
2471 f2(fTmp2),
2472 Ra1(Ra1Tmp),
2473 Ra2(Ra2Tmp),
2474 Ra1_3(Ra1Tmp.GetVec(3)),
2475 Ra2_3(Ra2Tmp.GetVec(3))
2476 {
2477 DEBUGCOUTFNAME("AerodynamicBeam2::AerodynamicBeam2");
2478
2479 ASSERT(pBeam != 0);
2480 ASSERT(pBeam->GetElemType() == Elem::BEAM);
2481
2482 pNode[NODE1] = pBeam->pGetNode(1);
2483 pNode[NODE2] = pBeam->pGetNode(2);
2484
2485 ASSERT(pNode[NODE1] != 0);
2486 ASSERT(pNode[NODE1]->GetNodeType() == Node::STRUCTURAL);
2487 ASSERT(pNode[NODE2] != 0);
2488 ASSERT(pNode[NODE2]->GetNodeType() == Node::STRUCTURAL);
2489 }
2490
~AerodynamicBeam2(void)2491 AerodynamicBeam2::~AerodynamicBeam2(void)
2492 {
2493 DEBUGCOUTFNAME("AerodynamicBeam2::~AerodynamicBeam2");
2494 }
2495
2496 /* Scrive il contributo dell'elemento al file di restart */
2497 std::ostream&
Restart(std::ostream & out) const2498 AerodynamicBeam2::Restart(std::ostream& out) const
2499 {
2500 DEBUGCOUTFNAME("AerodynamicBeam2::Restart");
2501 out << " aerodynamic beam2: " << GetLabel()
2502 << ", " << pBeam->GetLabel();
2503 if (pIndVel != 0) {
2504 out << ", rotor, " << pIndVel->GetLabel();
2505 }
2506 out << ", reference, node, ", f1.Write(out, ", ")
2507 << ", reference, node, 1, ", (Ra1.GetVec(1)).Write(out, ", ")
2508 << ", 2, ", (Ra1.GetVec(2)).Write(out, ", ")
2509 << ", reference, node, ", f2.Write(out, ", ")
2510 << ", reference, node, 1, ", (Ra2.GetVec(1)).Write(out, ", ")
2511 << ", 2, ", (Ra2.GetVec(2)).Write(out, ", ")
2512 << ", ";
2513 Chord.pGetShape()->Restart(out) << ", ";
2514 ForcePoint.pGetShape()->Restart(out) << ", ";
2515 VelocityPoint.pGetShape()->Restart(out) << ", ";
2516 Twist.pGetShape()->Restart(out) << ", "
2517 << GDI.iGetNum() << ", control, ";
2518 pGetDriveCaller()->Restart(out) << ", ";
2519 aerodata->Restart(out);
2520 return out << ";" << std::endl;
2521 }
2522
2523 static const doublereal pdsi2[] = { -1., 0. };
2524 static const doublereal pdsf2[] = { 0., 1. };
2525
2526 /* Jacobian assembly */
2527 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)2528 AerodynamicBeam2::AssJac(VariableSubMatrixHandler& WorkMat,
2529 doublereal dCoef,
2530 const VectorHandler& XCurr,
2531 const VectorHandler& XPrimeCurr)
2532 {
2533 DEBUGCOUT("Entering AerodynamicBeam2::AssJac()" << std::endl);
2534
2535 if (!bJacobian) {
2536 WorkMat.SetNullMatrix();
2537 return WorkMat;
2538 }
2539
2540 FullSubMatrixHandler& WM = WorkMat.SetFull();
2541
2542 /* Ridimensiona la sottomatrice in base alle esigenze */
2543 integer iNumRows = 0;
2544 integer iNumCols = 0;
2545 WorkSpaceDim(&iNumRows, &iNumCols);
2546 WM.ResizeReset(iNumRows, iNumCols);
2547
2548 integer iNode1FirstIndex = pNode[NODE1]->iGetFirstMomentumIndex();
2549 integer iNode1FirstPosIndex = pNode[NODE1]->iGetFirstPositionIndex();
2550 integer iNode2FirstIndex = pNode[NODE2]->iGetFirstMomentumIndex();
2551 integer iNode2FirstPosIndex = pNode[NODE2]->iGetFirstPositionIndex();
2552
2553 for (int iCnt = 1; iCnt <= 6; iCnt++) {
2554 WM.PutRowIndex(iCnt, iNode1FirstIndex + iCnt);
2555 WM.PutColIndex(iCnt, iNode1FirstPosIndex + iCnt);
2556 WM.PutRowIndex(6 + iCnt, iNode2FirstIndex + iCnt);
2557 WM.PutColIndex(6 + iCnt, iNode2FirstPosIndex + iCnt);
2558 }
2559
2560 doublereal dW[6];
2561
2562 /* array di vettori per via del ciclo sui nodi ... */
2563 Vec3 Xn[2];
2564
2565 /* Dati dei nodi */
2566 Xn[NODE1] = pNode[NODE1]->GetXCurr();
2567 const Mat3x3& Rn1(pNode[NODE1]->GetRCurr());
2568 const Vec3& Vn1(pNode[NODE1]->GetVCurr());
2569 const Vec3& Wn1(pNode[NODE1]->GetWCurr());
2570
2571 Xn[NODE2] = pNode[NODE2]->GetXCurr();
2572 const Mat3x3& Rn2(pNode[NODE2]->GetRCurr());
2573 const Vec3& Vn2(pNode[NODE2]->GetVCurr());
2574 const Vec3& Wn2(pNode[NODE2]->GetWCurr());
2575
2576 Vec3 f1Tmp(Rn1*f1);
2577 Vec3 f2Tmp(Rn2*f2);
2578
2579 Vec3 X1Tmp(Xn[NODE1] + f1Tmp);
2580 Vec3 X2Tmp(Xn[NODE2] + f2Tmp);
2581
2582 Vec3 V1Tmp(Vn1 + Wn1.Cross(f1Tmp));
2583 Vec3 V2Tmp(Vn2 + Wn2.Cross(f2Tmp));
2584
2585 Vec3 Omega1Crossf1(Wn1.Cross(f1Tmp));
2586 Vec3 Omega2Crossf2(Wn2.Cross(f2Tmp));
2587
2588 /*
2589 * Matrice di trasformazione dal sistema globale a quello aerodinamico
2590 */
2591 Mat3x3 RR1(Rn1*Ra1);
2592 Mat3x3 RR2(Rn2*Ra2);
2593
2594 /*
2595 * half of relative rotation vector from node 1 to node 2
2596 */
2597 Vec3 overline_theta(Vec3(ER_Rot::Param, RR1.MulTM(RR2))/2.);
2598
2599 /*
2600 * Se l'elemento e' collegato ad un rotore,
2601 * si fa dare la velocita' di rotazione
2602 */
2603 doublereal dOmega = 0.;
2604 if (pIndVel != 0) {
2605 Rotor *pRotor = dynamic_cast<Rotor *>(pIndVel);
2606 if (pRotor != 0) {
2607 dOmega = pRotor->dGetOmega();
2608 }
2609 }
2610
2611 /*
2612 * Dati "permanenti" (uso solo la posizione del nodo 2 perche'
2613 * non dovrebbero cambiare "molto")
2614 */
2615 Vec3 Xmid = (Xn[NODE2] + Xn[NODE1])/2.;
2616 doublereal rho, c, p, T;
2617 GetAirProps(Xmid, rho, c, p, T); /* p, T no used yet */
2618 aerodata->SetAirData(rho, c);
2619
2620 int iPnt = 0;
2621
2622 ResetIterator();
2623
2624 integer iNumDof = aerodata->iGetNumDof();
2625 integer iFirstEq = -1;
2626 integer iFirstSubEq = -1;
2627 if (iNumDof > 0) {
2628 iFirstEq = iGetFirstIndex();
2629 iFirstSubEq = 12;
2630
2631 integer iOffset = iFirstEq - 12;
2632
2633 for (int iCnt = 12 + 1; iCnt <= iNumRows; iCnt++) {
2634 WM.PutRowIndex(iCnt, iOffset + iCnt);
2635 WM.PutColIndex(iCnt, iOffset + iCnt);
2636 }
2637 }
2638
2639 for (int iNode = 0; iNode < LASTNODE; iNode++) {
2640 doublereal dsi = pdsi2[iNode];
2641 doublereal dsf = pdsf2[iNode];
2642
2643 doublereal dsm = (dsf + dsi)/2.;
2644 doublereal dsdCsi = (dsf - dsi)/2.;
2645
2646 Mat3x3 WM_F[4], WM_M[4];
2647 WM_F[DELTAx1].Reset();
2648 WM_F[DELTAg1].Reset();
2649 WM_F[DELTAx2].Reset();
2650 WM_F[DELTAg2].Reset();
2651 WM_M[DELTAx1].Reset();
2652 WM_M[DELTAg1].Reset();
2653 WM_M[DELTAx2].Reset();
2654 WM_M[DELTAg2].Reset();
2655
2656 /* Ciclo sui punti di Gauss */
2657 PntWght PW = GDI.GetFirst();
2658 do {
2659 doublereal dCsi = PW.dGetPnt();
2660 doublereal ds = dsm + dsdCsi*dCsi;
2661 doublereal dXds = DxDcsi2N(ds,
2662 Xn[NODE1], Xn[NODE2]);
2663
2664 doublereal dN1 = ShapeFunc2N(ds, 1);
2665 doublereal dN2 = ShapeFunc2N(ds, 2);
2666
2667 // note: identical to dN1 and dN2; see tecman.pdf
2668 #if 0
2669 doublereal dNN1 = (1. + dN1 - dN2)/2.;
2670 doublereal dNN2 = (1. + dN2 - dN1)/2.;
2671 #endif
2672
2673 Vec3 Xr(X1Tmp*dN1 + X2Tmp*dN2);
2674 Vec3 Vr(V1Tmp*dN1 + V2Tmp*dN2);
2675 Vec3 Wr(Wn1*dN1 + Wn2*dN2);
2676 Vec3 thetar(overline_theta*dN2);
2677
2678 /* Contributo di velocita' del vento */
2679 Vec3 VTmp(Zero3);
2680 if (fGetAirVelocity(VTmp, Xr)) {
2681 Vr -= VTmp;
2682 }
2683
2684 /*
2685 * Se l'elemento e' collegato ad un rotore,
2686 * aggiunge alla velocita' la velocita' indotta
2687 */
2688 if (pIndVel != 0) {
2689 Vr += pIndVel->GetInducedVelocity(GetElemType(),
2690 GetLabel(), iPnt, Xr);
2691 }
2692
2693 /* Copia i dati nel vettore di lavoro dVAM */
2694 doublereal dTw = Twist.dGet(ds);
2695 /* Contributo dell'eventuale sup. mobile */
2696 dTw += dGet();
2697
2698 aerodata->SetSectionData(dCsi,
2699 Chord.dGet(ds),
2700 ForcePoint.dGet(ds),
2701 VelocityPoint.dGet(ds),
2702 dTw,
2703 dOmega);
2704
2705 /*
2706 * Lo svergolamento non viene piu' trattato in aerod2_;
2707 * quindi lo uso per correggere la matrice di rotazione
2708 * dal sistema aerodinamico a quello globale
2709 */
2710 Mat3x3 RRloc(RR1*Mat3x3(ER_Rot::MatR, thetar));
2711 if (dTw != 0.) {
2712 doublereal dCosT = cos(dTw);
2713 doublereal dSinT = sin(dTw);
2714 /* Assumo lo svergolamento positivo a cabrare */
2715 Mat3x3 RTw(dCosT, dSinT, 0.,
2716 -dSinT, dCosT, 0.,
2717 0., 0., 1.);
2718 /*
2719 * Allo stesso tempo interpola le g
2720 * e aggiunge lo svergolamento
2721 */
2722 RRloc = RRloc*RTw;
2723 }
2724
2725 /*
2726 * Ruota velocita' e velocita' angolare nel sistema
2727 * aerodinamico e li copia nel vettore di lavoro dW
2728 */
2729 VTmp = RRloc.MulTV(Vr);
2730 VTmp.PutTo(&dW[0]);
2731
2732 Vec3 WTmp = RRloc.MulTV(Wr);
2733 WTmp.PutTo(&dW[3]);
2734 /* Funzione di calcolo delle forze aerodinamiche */
2735 doublereal Fa0[6];
2736 Mat6x6 JFa;
2737
2738 doublereal cc = dXds*dsdCsi*PW.dGetWght();
2739
2740 Vec3 d(Xr - Xn[iNode]);
2741
2742 Mat3x3 Bv1(MatCross, (Vr - Omega1Crossf1)*(dN1*dCoef));
2743 Mat3x3 Bv2(MatCross, (Vr - Omega2Crossf2)*(dN2*dCoef));
2744
2745 Mat3x3 Bw1(MatCross, (Wr - Wn1)*(dN1*dCoef));
2746 Mat3x3 Bw2(MatCross, (Wr - Wn2)*(dN2*dCoef));
2747
2748 if (iNumDof) {
2749 // prepare (v/dot{x} + dCoef*v/x) and so
2750 Mat3x3 RRlocT(RRloc.Transpose());
2751
2752 vx.PutMat3x3(1, RRlocT*dN1);
2753 vx.PutMat3x3(4, RRloc.MulTM(Bv1 - Mat3x3(MatCross, f1Tmp*dN1)));
2754
2755 vx.PutMat3x3(6 + 1, RRlocT*dN2);
2756 vx.PutMat3x3(6 + 4, RRloc.MulTM(Bv2 - Mat3x3(MatCross, f2Tmp*dN2)));
2757
2758 wx.PutMat3x3(4, RRlocT + Bw1);
2759 wx.PutMat3x3(6 + 4, RRlocT + Bw2);
2760
2761 // equations from iFirstEq on are dealt with by aerodata
2762 aerodata->AssJac(WM, dCoef, XCurr, XPrimeCurr,
2763 iFirstEq, iFirstSubEq,
2764 vx, wx, fq, cq, iPnt, dW, Fa0, JFa, OUTA[iPnt]);
2765
2766 // deal with (f/dot{q} + dCoef*f/q) and so
2767 integer iOffset = 12 + iPnt*iNumDof;
2768 for (integer iCol = 1; iCol <= iNumDof; iCol++) {
2769 Vec3 fqTmp((RRloc*fq.GetVec(iCol))*cc);
2770 Vec3 cqTmp(d.Cross(fqTmp) + (RRloc*cq.GetVec(iCol))*cc);
2771
2772 WM.Sub(6*iNode + 1, iOffset + iCol, fqTmp);
2773 WM.Sub(6*iNode + 4, iOffset + iCol, cqTmp);
2774 }
2775
2776 // first equation
2777 iFirstEq += iNumDof;
2778 iFirstSubEq += iNumDof;
2779
2780 } else {
2781 aerodata->GetForcesJac(iPnt, dW, Fa0, JFa, OUTA[iPnt]);
2782 }
2783
2784 // rotate force, couple and Jacobian matrix in absolute frame
2785 Mat6x6 JFaR = MultRMRt(JFa, RRloc, cc);
2786
2787 // force and moment about the node
2788 Vec3 fTmp(RRloc*(Vec3(&Fa0[0])*dCoef));
2789 Vec3 cTmp(RRloc*(Vec3(&Fa0[3])*dCoef) + d.Cross(fTmp));
2790
2791 Mat3x3 WM_F2[4];
2792
2793 // f <-> x
2794 WM_F2[DELTAx1] = JFaR.GetMat11()*dN1;
2795
2796 WM_F2[DELTAx2] = JFaR.GetMat11()*dN2;
2797
2798 doublereal delta;
2799
2800 // c <-> x
2801 delta = (iNode == NODE1) ? 1. : 0.;
2802 WM_M[DELTAx1] += JFaR.GetMat21()*dN1 - Mat3x3(MatCross, fTmp*(dN1 - delta));
2803
2804 delta = (iNode == NODE2) ? 1. : 0.;
2805 WM_M[DELTAx2] += JFaR.GetMat21()*dN2 - Mat3x3(MatCross, fTmp*(dN2 - delta));
2806
2807 // f <-> g
2808 WM_F2[DELTAg1] = (JFaR.GetMat12() - JFaR.GetMat11()*Mat3x3(MatCross, f1Tmp))*dN1;
2809 WM_F2[DELTAg1] += JFaR.GetMat11()*Bv1 + JFaR.GetMat12()*Bw1;
2810 WM_F2[DELTAg1] -= Mat3x3(MatCross, fTmp*dN1);
2811
2812 WM_F2[DELTAg2] = (JFaR.GetMat12() - JFaR.GetMat11()*Mat3x3(MatCross, f2Tmp))*dN2;
2813 WM_F2[DELTAg2] += JFaR.GetMat11()*Bv2 + JFaR.GetMat12()*Bw2;
2814 WM_F2[DELTAg2] -= Mat3x3(MatCross, fTmp*dN2);
2815
2816 // c <-> g
2817 WM_M[DELTAg1] += (JFaR.GetMat22() - JFaR.GetMat21()*Mat3x3(MatCross, f1Tmp))*dN1;
2818 WM_M[DELTAg1] += JFaR.GetMat21()*Bv1 + JFaR.GetMat22()*Bw1;
2819 WM_M[DELTAg1] -= Mat3x3(MatCross, cTmp*dN1);
2820 WM_M[DELTAg1] += Mat3x3(MatCrossCross, fTmp, f1Tmp*dN1);
2821
2822 WM_M[DELTAg2] += (JFaR.GetMat22() - JFaR.GetMat21()*Mat3x3(MatCross, f2Tmp))*dN2;
2823 WM_M[DELTAg2] += JFaR.GetMat21()*Bv2 + JFaR.GetMat22()*Bw2;
2824 WM_M[DELTAg2] -= Mat3x3(MatCross, cTmp*dN2);
2825 WM_M[DELTAg2] += Mat3x3(MatCrossCross, fTmp, f2Tmp*dN2);
2826
2827 for (int iCnt = 0; iCnt < 2*LASTNODE; iCnt++) {
2828 WM_F[iCnt] += WM_F2[iCnt];
2829 WM_M[iCnt] += d.Cross(WM_F2[iCnt]);
2830 }
2831
2832 iPnt++;
2833
2834 } while (GDI.fGetNext(PW));
2835
2836 for (int iCnt = 0; iCnt < 2*LASTNODE; iCnt++) {
2837 WM.Sub(6*iNode + 1, 3*iCnt + 1, WM_F[iCnt]);
2838 WM.Sub(6*iNode + 4, 3*iCnt + 1, WM_M[iCnt]);
2839 }
2840 }
2841
2842 return WorkMat;
2843 }
2844
2845 /* assemblaggio residuo */
2846 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)2847 AerodynamicBeam2::AssRes(
2848 SubVectorHandler& WorkVec,
2849 doublereal dCoef,
2850 const VectorHandler& XCurr,
2851 const VectorHandler& XPrimeCurr)
2852 {
2853 DEBUGCOUTFNAME("AerodynamicBeam2::AssRes");
2854
2855 integer iNumRows;
2856 integer iNumCols;
2857 WorkSpaceDim(&iNumRows, &iNumCols);
2858 WorkVec.ResizeReset(iNumRows);
2859
2860 integer iNode1FirstIndex = pNode[NODE1]->iGetFirstMomentumIndex();
2861 integer iNode2FirstIndex = pNode[NODE2]->iGetFirstMomentumIndex();
2862 for (int iCnt = 1; iCnt <= 6; iCnt++) {
2863 WorkVec.PutRowIndex(iCnt, iNode1FirstIndex + iCnt);
2864 WorkVec.PutRowIndex(6 + iCnt, iNode2FirstIndex + iCnt);
2865 }
2866
2867 AssVec(WorkVec, dCoef, XCurr, XPrimeCurr);
2868
2869 return WorkVec;
2870 }
2871
2872 /* assemblaggio iniziale residuo */
2873 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)2874 AerodynamicBeam2::InitialAssRes( SubVectorHandler& WorkVec,
2875 const VectorHandler& XCurr)
2876 {
2877 DEBUGCOUTFNAME("AerodynamicBeam2::InitialAssRes");
2878 WorkVec.ResizeReset(12);
2879
2880 integer iNode1FirstIndex = pNode[NODE1]->iGetFirstPositionIndex();
2881 integer iNode2FirstIndex = pNode[NODE2]->iGetFirstPositionIndex();
2882 for (int iCnt = 1; iCnt <= 6; iCnt++) {
2883 WorkVec.PutRowIndex(iCnt, iNode1FirstIndex + iCnt);
2884 WorkVec.PutRowIndex(6 + iCnt, iNode2FirstIndex + iCnt);
2885 }
2886
2887 AssVec(WorkVec, 1., XCurr, XCurr);
2888
2889 return WorkVec;
2890 }
2891
2892 /* assemblaggio residuo */
2893 void
AssVec(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)2894 AerodynamicBeam2::AssVec(SubVectorHandler& WorkVec,
2895 doublereal dCoef,
2896 const VectorHandler& XCurr,
2897 const VectorHandler& XPrimeCurr)
2898 {
2899 DEBUGCOUTFNAME("AerodynamicBeam2::AssVec");
2900
2901 doublereal dTng[6];
2902 doublereal dW[6];
2903
2904 Vec3 Xn[LASTNODE];
2905
2906 /* Dati dei nodi */
2907 Xn[NODE1] = pNode[NODE1]->GetXCurr();
2908 const Mat3x3& Rn1(pNode[NODE1]->GetRCurr());
2909 const Vec3& Vn1(pNode[NODE1]->GetVCurr());
2910 const Vec3& Wn1(pNode[NODE1]->GetWCurr());
2911
2912 Xn[NODE2] = pNode[NODE2]->GetXCurr();
2913 const Mat3x3& Rn2(pNode[NODE2]->GetRCurr());
2914 const Vec3& Vn2(pNode[NODE2]->GetVCurr());
2915 const Vec3& Wn2(pNode[NODE2]->GetWCurr());
2916
2917 Vec3 f1Tmp(Rn1*f1);
2918 Vec3 f2Tmp(Rn2*f2);
2919
2920 Vec3 X1Tmp(Xn[NODE1] + f1Tmp);
2921 Vec3 X2Tmp(Xn[NODE2] + f2Tmp);
2922
2923 Vec3 V1Tmp(Vn1 + Wn1.Cross(f1Tmp));
2924 Vec3 V2Tmp(Vn2 + Wn2.Cross(f2Tmp));
2925
2926 /*
2927 * Matrice di trasformazione dal sistema globale a quello aerodinamico
2928 */
2929 Mat3x3 RR1(Rn1*Ra1);
2930 Mat3x3 RR2(Rn2*Ra2);
2931
2932 /*
2933 * half of relative rotation vector from node 1 to node 2
2934 */
2935 Vec3 overline_theta(Vec3(ER_Rot::Param, RR1.MulTM(RR2))/2.);
2936
2937 /*
2938 * Se l'elemento e' collegato ad un rotore,
2939 * si fa dare la velocita' di rotazione
2940 */
2941 doublereal dOmega = 0.;
2942 if (pIndVel != 0) {
2943 Rotor *pRotor = dynamic_cast<Rotor *>(pIndVel);
2944 if (pRotor != 0) {
2945 dOmega = pRotor->dGetOmega();
2946 }
2947 }
2948
2949 /*
2950 * Dati "permanenti" (uso solo la posizione di mezzo perche'
2951 * non dovrebbero cambiare "molto")
2952 */
2953 Vec3 Xmid = (Xn[NODE2] + Xn[NODE1])/2.;
2954 doublereal rho, c, p, T;
2955 GetAirProps(Xmid, rho, c, p, T); /* p, T no used yet */
2956 aerodata->SetAirData(rho, c);
2957
2958 int iPnt = 0;
2959
2960 ResetIterator();
2961
2962 integer iNumDof = aerodata->iGetNumDof();
2963 integer iFirstEq = -1;
2964 integer iFirstSubEq = -1;
2965 if (iNumDof > 0) {
2966 iFirstEq = iGetFirstIndex();
2967 iFirstSubEq = 12;
2968
2969 integer iOffset = iFirstEq - 12;
2970 integer iNumRows = WorkVec.iGetSize();
2971 for (int iCnt = 12 + 1; iCnt <= iNumRows; iCnt++) {
2972 WorkVec.PutRowIndex(iCnt, iOffset + iCnt);
2973 }
2974 }
2975
2976 for (int iNode = 0; iNode < LASTNODE; iNode++) {
2977
2978 /* Resetta i dati */
2979 F[iNode].Reset();
2980 M[iNode].Reset();
2981
2982 doublereal dsi = pdsi2[iNode];
2983 doublereal dsf = pdsf2[iNode];
2984
2985 doublereal dsm = (dsf + dsi)/2.;
2986 doublereal dsdCsi = (dsf - dsi)/2.;
2987
2988 /* Ciclo sui punti di Gauss */
2989 PntWght PW = GDI.GetFirst();
2990 do {
2991 doublereal dCsi = PW.dGetPnt();
2992 doublereal ds = dsm + dsdCsi*dCsi;
2993 doublereal dXds = DxDcsi2N(ds, Xn[NODE1], Xn[NODE2]);
2994
2995 doublereal dN1 = ShapeFunc2N(ds, 1);
2996 doublereal dN2 = ShapeFunc2N(ds, 2);
2997
2998 Vec3 Xr(X1Tmp*dN1 + X2Tmp*dN2);
2999 Vec3 Vr(V1Tmp*dN1 + V2Tmp*dN2);
3000 Vec3 Wr(Wn1*dN1 + Wn2*dN2);
3001 Vec3 thetar(overline_theta*((1. + dN2 - dN1)/2.));
3002
3003 /* Contributo di velocita' del vento */
3004 Vec3 VTmp(Zero3);
3005 if (fGetAirVelocity(VTmp, Xr)) {
3006 Vr -= VTmp;
3007 }
3008
3009 /*
3010 * Se l'elemento e' collegato ad un rotore,
3011 * aggiunge alla velocita' la velocita' indotta
3012 */
3013 if (pIndVel != 0) {
3014 Vr += pIndVel->GetInducedVelocity(GetElemType(),
3015 GetLabel(), iPnt, Xr);
3016 }
3017
3018 /* Copia i dati nel vettore di lavoro dVAM */
3019 doublereal dTw = Twist.dGet(ds);
3020 dTw += dGet(); /* Contributo dell'eventuale sup. mobile */
3021 aerodata->SetSectionData(dCsi,
3022 Chord.dGet(ds),
3023 ForcePoint.dGet(ds),
3024 VelocityPoint.dGet(ds),
3025 dTw,
3026 dOmega);
3027
3028 /*
3029 * Lo svergolamento non viene piu' trattato in aerod2_; quindi
3030 * lo uso per correggere la matrice di rotazione
3031 * dal sistema aerodinamico a quello globale
3032 */
3033 Mat3x3 RRloc(RR1*Mat3x3(ER_Rot::MatR, thetar));
3034 if (dTw != 0.) {
3035 doublereal dCosT = cos(dTw);
3036 doublereal dSinT = sin(dTw);
3037 /* Assumo lo svergolamento positivo a cabrare */
3038 Mat3x3 RTw( dCosT, dSinT, 0.,
3039 -dSinT, dCosT, 0.,
3040 0., 0., 1.);
3041 /*
3042 * Allo stesso tempo interpola le g e aggiunge lo svergolamento
3043 */
3044 RRloc = RRloc*RTw;
3045 }
3046
3047 /*
3048 * Ruota velocita' e velocita' angolare nel sistema
3049 * aerodinamico e li copia nel vettore di lavoro dW
3050 */
3051 VTmp = RRloc.MulTV(Vr);
3052 VTmp.PutTo(dW);
3053
3054 Vec3 WTmp = RRloc.MulTV(Wr);
3055 WTmp.PutTo(&dW[3]);
3056
3057 /* Funzione di calcolo delle forze aerodinamiche */
3058 if (iNumDof) {
3059 aerodata->AssRes(WorkVec, dCoef, XCurr, XPrimeCurr,
3060 iFirstEq, iFirstSubEq, iPnt, dW, dTng, OUTA[iPnt]);
3061
3062 // first equation
3063 iFirstEq += iNumDof;
3064 iFirstSubEq += iNumDof;
3065
3066 } else {
3067 aerodata->GetForces(iPnt, dW, dTng, OUTA[iPnt]);
3068 }
3069
3070 /* Dimensionalizza le forze */
3071 doublereal dWght = dXds*dsdCsi*PW.dGetWght();
3072 dTng[1] *= TipLoss.dGet(dCsi);
3073 Vec3 FTmp(RRloc*(Vec3(&dTng[0])));
3074 Vec3 MTmp(RRloc*(Vec3(&dTng[3])));
3075
3076 // Se e' definito il rotore, aggiungere il contributo alla trazione
3077 AddSectionalForce_int(iPnt, FTmp, MTmp, dWght, Xr, RRloc, Vr, Wr);
3078
3079 FTmp *= dWght;
3080 MTmp *= dWght;
3081 F[iNode] += FTmp;
3082 M[iNode] += MTmp;
3083 M[iNode] += (Xr - Xn[iNode]).Cross(FTmp);
3084
3085 // specific for Gauss points force output
3086 if (bToBeOutput()) {
3087 SetData(VTmp, dTng, Xr, RRloc, Vr, Wr, FTmp, MTmp);
3088 }
3089
3090 iPnt++;
3091
3092 } while (GDI.fGetNext(PW));
3093
3094 // Se e' definito il rotore, aggiungere il contributo alla trazione
3095 AddForce_int(pNode[iNode], F[iNode], M[iNode], Xn[iNode]);
3096
3097 /* Somma il termine al residuo */
3098 WorkVec.Add(6*iNode + 1, F[iNode]);
3099 WorkVec.Add(6*iNode + 4, M[iNode]);
3100 }
3101 }
3102
3103 /*
3104 * output; si assume che ogni tipo di elemento sappia, attraverso
3105 * l'OutputHandler, dove scrivere il proprio output
3106 */
3107 void
Output(OutputHandler & OH) const3108 AerodynamicBeam2::Output(OutputHandler& OH ) const
3109 {
3110 DEBUGCOUTFNAME("AerodynamicBeam2::Output");
3111
3112 if (bToBeOutput()) {
3113 Aerodynamic2DElem<2>::Output_int(OH);
3114
3115 if (OH.UseText(OutputHandler::AERODYNAMIC)) {
3116 std::ostream& out = OH.Aerodynamic() << std::setw(8) << GetLabel();
3117
3118 switch (GetOutput()) {
3119
3120 case AEROD_OUT_NODE:
3121 out << " " << std::setw(8) << pBeam->GetLabel()
3122 << " ", F[NODE1].Write(out, " ") << " ", M[NODE1].Write(out, " ")
3123 << " ", F[NODE2].Write(out, " ") << " ", M[NODE2].Write(out, " ");
3124 break;
3125
3126 case AEROD_OUT_PGAUSS:
3127 ASSERT(!OutputData.empty());
3128
3129 for (std::vector<Aero_output>::const_iterator i = OutputData.begin();
3130 i != OutputData.end(); ++i)
3131 {
3132 out << " " << i->alpha
3133 << " " << i->f;
3134 }
3135 break;
3136
3137 case AEROD_OUT_STD:
3138 for (int i = 0; i < 2*GDI.iGetNum(); i++) {
3139 out
3140 << " " << OUTA[i].alpha
3141 << " " << OUTA[i].gamma
3142 << " " << OUTA[i].mach
3143 << " " << OUTA[i].cl
3144 << " " << OUTA[i].cd
3145 << " " << OUTA[i].cm
3146 << " " << OUTA[i].alf1
3147 << " " << OUTA[i].alf2;
3148 }
3149 break;
3150
3151 default:
3152 ASSERT(0);
3153 break;
3154 }
3155
3156 out << std::endl;
3157 }
3158 }
3159 }
3160
3161 /* AerodynamicBeam2 - end */
3162
3163
3164 /* Legge un elemento aerodinamico di trave a due nodi */
3165
3166 Elem *
ReadAerodynamicBeam2(DataManager * pDM,MBDynParser & HP,const DofOwner * pDO,unsigned int uLabel)3167 ReadAerodynamicBeam2(DataManager* pDM,
3168 MBDynParser& HP,
3169 const DofOwner *pDO,
3170 unsigned int uLabel)
3171 {
3172 DEBUGCOUTFNAME("ReadAerodynamicBeam2");
3173
3174 /* Trave */
3175 unsigned int uBeam = (unsigned int)HP.GetInt();
3176
3177 DEBUGLCOUT(MYDEBUG_INPUT, "Linked to beam: " << uBeam << std::endl);
3178
3179 /* verifica di esistenza della trave */
3180 Elem *p = pDM->pFindElem(Elem::BEAM, uBeam);
3181 if (p == 0) {
3182 silent_cerr("Beam2(" << uBeam << ") not defined "
3183 "at line " << HP.GetLineData()
3184 << std::endl);
3185 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3186 }
3187 Beam2* pBeam = dynamic_cast<Beam2 *>(p);
3188 if (pBeam == 0) {
3189 silent_cerr("Beam(" << uBeam << ") is not a Beam2 "
3190 "at line " << HP.GetLineData()
3191 << std::endl);
3192 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3193 }
3194
3195 /* Eventuale rotore */
3196 InducedVelocity* pIndVel = 0;
3197 bool bPassive(false);
3198 (void)ReadInducedVelocity(pDM, HP, uLabel, "AerodynamicBeam2",
3199 pIndVel, bPassive);
3200
3201 /* Nodo 1: */
3202
3203 /* Offset del corpo aerodinamico rispetto al nodo */
3204 const StructNode* pNode1 = pBeam->pGetNode(1);
3205
3206 ReferenceFrame RF(pNode1);
3207 Vec3 f1(HP.GetPosRel(RF));
3208 DEBUGLCOUT(MYDEBUG_INPUT, "Node 1 offset: " << f1 << std::endl);
3209
3210 Mat3x3 Ra1(HP.GetRotRel(RF));
3211 DEBUGLCOUT(MYDEBUG_INPUT,
3212 "Node 1 rotation matrix: " << std::endl << Ra1 << std::endl);
3213
3214 /* Nodo 2: */
3215
3216 /* Offset del corpo aerodinamico rispetto al nodo */
3217 const StructNode* pNode2 = pBeam->pGetNode(2);
3218
3219 RF = ReferenceFrame(pNode2);
3220 Vec3 f2(HP.GetPosRel(RF));
3221 DEBUGLCOUT(MYDEBUG_INPUT, "Node 2 offset: " << f2 << std::endl);
3222
3223 Mat3x3 Ra2(HP.GetRotRel(RF));
3224 DEBUGLCOUT(MYDEBUG_INPUT,
3225 "Node 2 rotation matrix: " << std::endl << Ra2 << std::endl);
3226
3227 Shape* pChord = 0;
3228 Shape* pForce = 0;
3229 Shape* pVelocity = 0;
3230 Shape* pTwist = 0;
3231 Shape* pTipLoss = 0;
3232
3233 integer iNumber = 0;
3234 DriveCaller* pDC = 0;
3235 AeroData* aerodata = 0;
3236
3237 ReadAeroData(pDM, HP, 2,
3238 &pChord, &pForce, &pVelocity, &pTwist, &pTipLoss,
3239 &iNumber, &pDC, &aerodata);
3240
3241 bool bUseJacobian(false);
3242 if (HP.IsKeyWord("jacobian")) {
3243 bUseJacobian = HP.GetYesNoOrBool(bDefaultUseJacobian);
3244 }
3245
3246 if (aerodata->iGetNumDof() > 0 && !bUseJacobian) {
3247 silent_cerr("AerodynamicBeam2(" << uLabel << "): "
3248 "aerodynamic model needs \"jacobian, yes\" at line " << HP.GetLineData()
3249 << std::endl);
3250 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3251 }
3252
3253 OrientationDescription od = UNKNOWN_ORIENTATION_DESCRIPTION;
3254 unsigned uFlags = AerodynamicOutput::OUTPUT_NONE;
3255 ReadOptionalAerodynamicCustomOutput(pDM, HP, uLabel, uFlags, od);
3256
3257 flag fOut = pDM->fReadOutput(HP, Elem::AERODYNAMIC);
3258 if (HP.IsArg()) {
3259 if (HP.IsKeyWord("std")) {
3260 fOut |= AerodynamicOutput::AEROD_OUT_STD;
3261 } else if (HP.IsKeyWord("gauss")) {
3262 fOut |= AerodynamicOutput::AEROD_OUT_PGAUSS;
3263 } else if (HP.IsKeyWord("node")) {
3264 fOut |= AerodynamicOutput::AEROD_OUT_NODE;
3265 } else {
3266 silent_cerr("AerodynamicBeam2(" << uLabel << "): "
3267 "unknown output mode at line " << HP.GetLineData()
3268 << std::endl);
3269 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3270 }
3271
3272 } else if (fOut) {
3273 fOut |= AerodynamicOutput::AEROD_OUT_STD;
3274 }
3275 fOut |= uFlags;
3276
3277 Elem* pEl = 0;
3278
3279 SAFENEWWITHCONSTRUCTOR(pEl,
3280 AerodynamicBeam2,
3281 AerodynamicBeam2(uLabel, pDO, pBeam, pIndVel, bPassive,
3282 f1, f2, Ra1, Ra2,
3283 pChord, pForce, pVelocity, pTwist, pTipLoss,
3284 iNumber, aerodata, pDC, bUseJacobian, od, fOut));
3285
3286 /* Se non c'e' il punto e virgola finale */
3287 if (HP.IsArg()) {
3288 silent_cerr("semicolon expected at line "
3289 << HP.GetLineData() << std::endl);
3290 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3291 }
3292
3293 std::ostream& out = pDM->GetLogFile();
3294 out << "aero2: " << uLabel;
3295
3296 Vec3 ra1 = Ra1.GetVec(1);
3297 Vec3 ra3 = Ra1.GetVec(3);
3298 doublereal dC = pChord->dGet(-1.);
3299 doublereal dP = pForce->dGet(-1.);
3300 out
3301 << " " << pNode1->GetLabel()
3302 << " ", (f1 + ra1*(dP - dC*3./4.)).Write(out, " ")
3303 << " ", (f1 + ra1*(dP + dC/4.)).Write(out, " ");
3304
3305 ra1 = Ra2.GetVec(1);
3306 ra3 = Ra2.GetVec(3);
3307 dC = pChord->dGet(1.);
3308 dP = pForce->dGet(1.);
3309 out
3310 << " " << pNode2->GetLabel()
3311 << " ", (f2 + ra1*(dP - dC*3./4.)).Write(out, " ")
3312 << " ", (f2 + ra1*(dP + dC/4.)).Write(out, " ")
3313 << std::endl;
3314
3315 return pEl;
3316 } /* End of ReadAerodynamicBeam2() */
3317
3318