1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/genfm.cc,v 1.33 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 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33
34 #include <cmath>
35
36 #include "dataman.h"
37 #include "genfm.h"
38 #include "bisec.h"
39
40 /* GenericAerodynamicForce - begin */
41
42 // TODO: add private data; e.g. alpha, beta
43
44 static const doublereal dAlphaMax[] = { M_PI/2, M_PI };
45 static const doublereal dBetaMax[] = { M_PI, M_PI/2 };
46
47 /* Assemblaggio residuo */
48 void
AssVec(SubVectorHandler & WorkVec)49 GenericAerodynamicForce::AssVec(SubVectorHandler& WorkVec)
50 {
51 /* velocity at aerodynamic center */
52
53 /* position of aerodynamic center */
54 const Mat3x3& Rn(pNode->GetRCurr());
55 Vec3 f(Rn*tilde_f);
56 Mat3x3 R(Rn*tilde_Ra);
57 Vec3 Xca(pNode->GetXCurr() + f);
58 Vec3 Vca(pNode->GetVCurr() + f.Cross(pNode->GetWCurr()));
59
60 Vec3 VTmp(Zero3);
61 if (fGetAirVelocity(VTmp, Xca)) {
62 Vca -= VTmp;
63 }
64
65 #if 0
66 /*
67 * Se l'elemento e' collegato ad un rotore,
68 * aggiunge alla velocita' la velocita' indotta
69 */
70 if (pRotor != NULL) {
71 Vca += pRotor->GetInducedVelocity(GetElemType(), GetLabel(), 0, Xca);
72 }
73 #endif
74
75 Vec3 V(R.MulTV(Vca));
76
77 if (bAlphaFirst) {
78 dAlpha = atan2(V(3), copysign(std::sqrt(V(1)*V(1) + V(2)*V(2)), V(1)));
79 dBeta = -atan2(V(2), std::abs(V(1)));
80
81 } else {
82 dAlpha = atan2(V(3), std::abs(V(1)));
83 dBeta = -atan2(V(2), copysign(std::sqrt(V(1)*V(1) + V(2)*V(2)), V(1)));
84 }
85
86 doublereal rho, c, p, T;
87 GetAirProps(Xca, rho, c, p, T); /* p, T no used yet */
88 doublereal q = Vca.Dot()*rho/2.;
89
90 doublereal dScaleForce = q*dRefSurface;
91 doublereal dScaleMoment = dScaleForce*dRefLength;
92
93 int nAlpha = pData->Alpha.size() - 1;
94 int nBeta = pData->Beta.size() - 1;
95
96 ASSERT(nAlpha >= 0);
97 ASSERT(nBeta >= 0);
98
99 if (dAlpha < pData->Alpha[0]) {
100 /* smooth out coefficients if Alpha does not span -180 => 180 */
101 doublereal dAlphaX = (dAlpha - pData->Alpha[0])/(-::dAlphaMax[bAlphaFirst] - pData->Alpha[0]);
102 doublereal dSmoothAlpha = (std::cos(::dAlphaMax[bAlphaFirst]*dAlphaX) + 1)/2.;
103
104 if (dBeta < pData->Beta[0]) {
105 /* smooth out coefficients if Beta does not span -180 => 180 */
106 doublereal dBetaX = (dBeta - pData->Beta[0])/(-::dBetaMax[bAlphaFirst] - pData->Beta[0]);
107 doublereal dSmoothBeta = (std::cos(::dBetaMax[bAlphaFirst]*dBetaX) + 1)/2.;
108
109 tilde_F = Vec3(&pData->Data[0][0].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
110 tilde_M = Vec3(&pData->Data[0][0].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
111
112 } else if (dBeta > pData->Beta[nBeta]) {
113 /* smooth out coefficients if Beta does not span -180 => 180 */
114 doublereal dBetaX = (dBeta - pData->Beta[nBeta])/(::dBetaMax[bAlphaFirst] - pData->Beta[nBeta]);
115 doublereal dSmoothBeta = (std::cos(::dBetaMax[bAlphaFirst]*dBetaX) + 1)/2.;
116
117 tilde_F = Vec3(&pData->Data[nBeta][0].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
118 tilde_M = Vec3(&pData->Data[nBeta][0].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
119
120 } else {
121 int iBeta = bisec<doublereal>(&pData->Beta[0], dBeta, 0, nBeta);
122
123 ASSERT(iBeta >= 0);
124 ASSERT(iBeta < nBeta);
125
126 doublereal ddBeta = pData->Beta[iBeta + 1] - pData->Beta[iBeta];
127 doublereal d1Beta = (pData->Beta[iBeta + 1] - dBeta)/ddBeta;
128 doublereal d2Beta = (dBeta - pData->Beta[iBeta])/ddBeta;
129
130 GenericAerodynamicData::GenericAerodynamicCoef c
131 = pData->Data[iBeta][0]*d1Beta + pData->Data[iBeta + 1][0]*d2Beta;
132
133 tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothAlpha);
134 tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothAlpha);
135 }
136
137 } else if (dAlpha > pData->Alpha[nAlpha]) {
138 /* smooth out coefficients if Alpha does not span -180 => 180 */
139 doublereal dAlphaX = (dAlpha - pData->Alpha[nAlpha])/(-::dAlphaMax[bAlphaFirst] - pData->Alpha[nAlpha]);
140 doublereal dSmoothAlpha = (std::cos(::dAlphaMax[bAlphaFirst]*dAlphaX) + 1)/2.;
141
142 if (dBeta < pData->Beta[0]) {
143 /* smooth out coefficients if Beta does not span -180 => 180 */
144 doublereal dBetaX = (dBeta - pData->Beta[0])/(-::dBetaMax[bAlphaFirst] - pData->Beta[0]);
145 doublereal dSmoothBeta = (std::cos(::dBetaMax[bAlphaFirst]*dBetaX) + 1)/2.;
146
147 tilde_F = Vec3(&pData->Data[0][nAlpha].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
148 tilde_M = Vec3(&pData->Data[0][nAlpha].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
149
150 } else if (dBeta > pData->Beta[nBeta]) {
151 /* smooth out coefficients if Beta does not span -180 => 180 */
152 doublereal dBetaX = (dBeta - pData->Beta[nBeta])/(::dBetaMax[bAlphaFirst] - pData->Beta[nBeta]);
153 doublereal dSmoothBeta = (std::cos(::dBetaMax[bAlphaFirst]*dBetaX) + 1)/2.;
154
155 tilde_F = Vec3(&pData->Data[nBeta][nAlpha].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
156 tilde_M = Vec3(&pData->Data[nBeta][nAlpha].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
157
158 } else {
159 int iBeta = bisec<doublereal>(&pData->Beta[0], dBeta, 0, nBeta);
160
161 ASSERT(iBeta >= 0);
162 ASSERT(iBeta < nBeta);
163
164 doublereal ddBeta = pData->Beta[iBeta + 1] - pData->Beta[iBeta];
165 doublereal d1Beta = (pData->Beta[iBeta + 1] - dBeta)/ddBeta;
166 doublereal d2Beta = (dBeta - pData->Beta[iBeta])/ddBeta;
167
168 GenericAerodynamicData::GenericAerodynamicCoef c
169 = pData->Data[iBeta][nAlpha]*d1Beta + pData->Data[iBeta + 1][nAlpha]*d2Beta;
170
171 tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothAlpha);
172 tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothAlpha);
173 }
174
175 } else {
176 int iAlpha = bisec<doublereal>(&pData->Alpha[0], dAlpha, 0, nAlpha);
177
178 ASSERT(iAlpha >= 0);
179 ASSERT(iAlpha < nAlpha);
180
181 doublereal ddAlpha = pData->Alpha[iAlpha + 1] - pData->Alpha[iAlpha];
182 doublereal d1Alpha = (pData->Alpha[iAlpha + 1] - dAlpha)/ddAlpha;
183 doublereal d2Alpha = (dAlpha - pData->Alpha[iAlpha])/ddAlpha;
184
185 if (dBeta < pData->Beta[0]) {
186 /* smooth out coefficients if Beta does not span -180 => 180 */
187 doublereal dBetaX = (dBeta - pData->Beta[0])/(-::dBetaMax[bAlphaFirst] - pData->Beta[0]);
188 doublereal dSmoothBeta = (std::cos(::dBetaMax[bAlphaFirst]*dBetaX) + 1)/2.;
189
190 GenericAerodynamicData::GenericAerodynamicCoef c
191 = pData->Data[0][iAlpha]*d1Alpha + pData->Data[0][iAlpha + 1]*d2Alpha;
192
193 tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothBeta);
194 tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothBeta);
195
196 } else if (dBeta > pData->Beta[nBeta]) {
197 /* smooth out coefficients if Beta does not span -180 => 180 */
198 doublereal dBetaX = (dBeta - pData->Beta[nBeta])/(::dBetaMax[bAlphaFirst] - pData->Beta[nBeta]);
199 doublereal dSmoothBeta = (std::cos(::dBetaMax[bAlphaFirst]*dBetaX) + 1)/2.;
200
201 GenericAerodynamicData::GenericAerodynamicCoef c
202 = pData->Data[nBeta][iAlpha]*d1Alpha + pData->Data[nBeta][iAlpha + 1]*d2Alpha;
203
204 tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothBeta);
205 tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothBeta);
206
207 } else {
208 int iBeta = bisec<doublereal>(&pData->Beta[0], dBeta, 0, nBeta);
209
210 ASSERT(iBeta >= 0);
211 ASSERT(iBeta < nBeta);
212
213 doublereal ddBeta = pData->Beta[iBeta + 1] - pData->Beta[iBeta];
214 doublereal d1Beta = (pData->Beta[iBeta + 1] - dBeta)/ddBeta;
215 doublereal d2Beta = (dBeta - pData->Beta[iBeta])/ddBeta;
216
217 GenericAerodynamicData::GenericAerodynamicCoef c1
218 = pData->Data[iBeta][iAlpha]*d1Alpha + pData->Data[iBeta][iAlpha + 1]*d2Alpha;
219 GenericAerodynamicData::GenericAerodynamicCoef c2
220 = pData->Data[iBeta + 1][iAlpha]*d1Alpha + pData->Data[iBeta + 1][iAlpha + 1]*d2Alpha;
221
222 GenericAerodynamicData::GenericAerodynamicCoef c = c1*d1Beta + c2*d2Beta;
223
224 tilde_F = Vec3(&c.dCoef[0])*dScaleForce;
225 tilde_M = Vec3(&c.dCoef[3])*dScaleMoment;
226 }
227 }
228
229 F = R*tilde_F;
230 M = R*tilde_M + f.Cross(F);
231
232 WorkVec.Add(1, F);
233 WorkVec.Add(4, M);
234 }
235
GenericAerodynamicForce(unsigned int uLabel,const DofOwner * pDO,const StructNode * pN,const Vec3 & fTmp,const Mat3x3 & RaTmp,doublereal dS,doublereal dL,bool bAlphaFirst,GenericAerodynamicData * pD,flag fOut)236 GenericAerodynamicForce::GenericAerodynamicForce(unsigned int uLabel,
237 const DofOwner *pDO,
238 const StructNode* pN,
239 const Vec3& fTmp, const Mat3x3& RaTmp,
240 doublereal dS, doublereal dL, bool bAlphaFirst,
241 GenericAerodynamicData *pD,
242 flag fOut)
243 : Elem(uLabel, fOut),
244 AerodynamicElem(uLabel, pDO, fOut),
245 InitialAssemblyElem(uLabel, fOut),
246 pNode(pN),
247 dRefSurface(dS),
248 dRefLength(dL),
249 bAlphaFirst(bAlphaFirst),
250 tilde_f(fTmp),
251 tilde_Ra(RaTmp),
252 tilde_F(Zero3),
253 tilde_M(Zero3),
254 F(Zero3),
255 M(Zero3),
256 dAlpha(0.),
257 dBeta(0.),
258 pData(pD)
259 {
260 NO_OP;
261 }
262
~GenericAerodynamicForce(void)263 GenericAerodynamicForce::~GenericAerodynamicForce(void)
264 {
265 if (pData != 0) {
266 delete pData;
267 }
268 }
269
270 /* Tipo dell'elemento (usato per debug ecc.) */
271 Elem::Type
GetElemType(void) const272 GenericAerodynamicForce::GetElemType(void) const
273 {
274 return Elem::AERODYNAMIC;
275 }
276
277 /* Dimensioni del workspace */
278 void
WorkSpaceDim(integer * piNumRows,integer * piNumCols) const279 GenericAerodynamicForce::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
280 {
281 *piNumRows = 6;
282 *piNumCols = 1;
283 }
284
285 /* assemblaggio jacobiano */
286 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal,const VectorHandler &,const VectorHandler &)287 GenericAerodynamicForce::AssJac(VariableSubMatrixHandler& WorkMat,
288 doublereal /* dCoef */ ,
289 const VectorHandler& /* XCurr */ ,
290 const VectorHandler& /* XPrimeCurr */ )
291 {
292 DEBUGCOUTFNAME("GenericAerodynamicForce::AssJac");
293 WorkMat.SetNullMatrix();
294 return WorkMat;
295 }
296
297 /* Numero di GDL iniziali */
298 unsigned int
iGetInitialNumDof(void) const299 GenericAerodynamicForce::iGetInitialNumDof(void) const
300 {
301 return 0;
302 }
303
304 /* Dimensioni del workspace */
305 void
InitialWorkSpaceDim(integer * piNumRows,integer * piNumCols) const306 GenericAerodynamicForce::InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const
307 {
308 *piNumRows = 6;
309 *piNumCols = 1;
310 }
311
312 /* assemblaggio jacobiano */
313 VariableSubMatrixHandler&
InitialAssJac(VariableSubMatrixHandler & WorkMat,const VectorHandler &)314 GenericAerodynamicForce::InitialAssJac(VariableSubMatrixHandler& WorkMat,
315 const VectorHandler& /* XCurr */)
316 {
317 DEBUGCOUTFNAME("GenericAerodynamicForce::InitialAssJac");
318 WorkMat.SetNullMatrix();
319 return WorkMat;
320 }
321
322 /* Tipo di elemento aerodinamico */
323 AerodynamicElem::Type
GetAerodynamicElemType(void) const324 GenericAerodynamicForce::GetAerodynamicElemType(void) const
325 {
326 return AerodynamicElem::GENERICFORCE;
327 }
328
329 void
GetConnectedNodes(std::vector<const Node * > & connectedNodes) const330 GenericAerodynamicForce::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
331 {
332 connectedNodes.resize(1);
333 connectedNodes[0] = pNode;
334 }
335
336 /* Scrive il contributo dell'elemento al file di restart */
337 std::ostream&
Restart(std::ostream & out) const338 GenericAerodynamicForce::Restart(std::ostream& out) const
339 {
340 return out;
341 }
342
343 /* assemblaggio residuo */
344 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr)345 GenericAerodynamicForce::AssRes(SubVectorHandler& WorkVec,
346 doublereal dCoef,
347 const VectorHandler& XCurr,
348 const VectorHandler& XPrimeCurr)
349 {
350 WorkVec.ResizeReset(6);
351
352 integer iNodeFirstIndex = pNode->iGetFirstMomentumIndex();
353 for (int iCnt = 1; iCnt <= 6; iCnt++) {
354 WorkVec.PutRowIndex(iCnt, iNodeFirstIndex + iCnt);
355 }
356
357 AssVec(WorkVec);
358
359 return WorkVec;
360 }
361
362 /*
363 * output; si assume che ogni tipo di elemento sappia, attraverso
364 * l'OutputHandler, dove scrivere il proprio output
365 */
366 void
Output(OutputHandler & OH) const367 GenericAerodynamicForce::Output(OutputHandler& OH) const
368 {
369 if (bToBeOutput()) {
370 OH.Aerodynamic()
371 << std::setw(8) << GetLabel()
372 << " " << dAlpha*180./M_PI
373 << " " << dBeta*180./M_PI
374 << " " << tilde_F << " " << tilde_M
375 << " " << F << " " << M << std::endl;
376 }
377 }
378
379 /* Dati privati */
380 unsigned int
iGetNumPrivData(void) const381 GenericAerodynamicForce::iGetNumPrivData(void) const
382 {
383 return 2 + 6 + 6;
384 }
385
386 unsigned int
iGetPrivDataIdx(const char * s) const387 GenericAerodynamicForce::iGetPrivDataIdx(const char *s) const
388 {
389 unsigned idx = 0;
390 switch (s[0]) {
391 case 'F':
392 break;
393
394 case 'M':
395 idx += 3;
396 break;
397
398 case 'f':
399 idx += 6;
400 break;
401
402 case 'm':
403 idx += 9;
404 break;
405
406 default:
407 if (strcmp(s, "alpha") == 0) {
408 return idx + 12 + 1;
409
410 } else if (strcmp(s, "beta") == 0) {
411 return idx + 12 + 2;
412
413 }
414 return 0;
415 }
416
417 if (s[2] != '\0') {
418 return 0;
419 }
420
421 switch (s[1]) {
422 case 'x':
423 idx += 1;
424 break;
425
426 case 'y':
427 idx += 2;
428 break;
429
430 case 'z':
431 idx += 3;
432 break;
433
434 default:
435 return 0;
436 }
437
438 return idx;
439 }
440
441 doublereal
dGetPrivData(unsigned int i) const442 GenericAerodynamicForce::dGetPrivData(unsigned int i) const
443 {
444 if (i <= 0 || i > iGetNumPrivData()) {
445 // error
446 return 0.;
447 }
448
449 switch (i) {
450 case 1:
451 case 2:
452 case 3:
453 return F(i);
454
455 case 3 + 1:
456 case 3 + 2:
457 case 3 + 3:
458 return M(i - 3);
459
460 case 6 + 1:
461 case 6 + 2:
462 case 6 + 3:
463 return tilde_F(i - 6);
464
465 case 9 + 1:
466 case 9 + 2:
467 case 9 + 3:
468 return tilde_M(i - 9);
469
470 case 12 + 1:
471 return dAlpha;
472
473 case 12 + 2:
474 return dBeta;
475
476 default:
477 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
478 }
479 }
480
481 /* assemblaggio residuo */
482 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)483 GenericAerodynamicForce::InitialAssRes(SubVectorHandler& WorkVec,
484 const VectorHandler& XCurr)
485 {
486 return WorkVec;
487 }
488
489 extern Elem *
490 ReadGenericAerodynamicForce(DataManager* pDM,
491 MBDynParser& HP,
492 unsigned int uLabel);
493
494 /* GenericAerodynamicForce - end */
495
496 /* GenericAerodynamicData - begin */
497
498 /* GenericAerodynamicData::GenericAerodynamicCoef - begin */
499
GenericAerodynamicCoef(void)500 GenericAerodynamicData::GenericAerodynamicCoef::GenericAerodynamicCoef(void)
501 {
502 memset(&dCoef[0], 0, sizeof(dCoef));
503 }
504
GenericAerodynamicCoef(const GenericAerodynamicData::GenericAerodynamicCoef & c)505 GenericAerodynamicData::GenericAerodynamicCoef::GenericAerodynamicCoef(
506 const GenericAerodynamicData::GenericAerodynamicCoef& c)
507 {
508 dCoef[0] = c.dCoef[0];
509 dCoef[1] = c.dCoef[1];
510 dCoef[2] = c.dCoef[2];
511 dCoef[3] = c.dCoef[3];
512 dCoef[4] = c.dCoef[4];
513 dCoef[5] = c.dCoef[5];
514 }
515
516 GenericAerodynamicData::GenericAerodynamicCoef
operator +(const GenericAerodynamicData::GenericAerodynamicCoef & c) const517 GenericAerodynamicData::GenericAerodynamicCoef::operator + (
518 const GenericAerodynamicData::GenericAerodynamicCoef& c) const
519 {
520 GenericAerodynamicData::GenericAerodynamicCoef retval;
521
522 retval.dCoef[0] = dCoef[0] + c.dCoef[0];
523 retval.dCoef[1] = dCoef[1] + c.dCoef[1];
524 retval.dCoef[2] = dCoef[2] + c.dCoef[2];
525 retval.dCoef[3] = dCoef[3] + c.dCoef[3];
526 retval.dCoef[4] = dCoef[4] + c.dCoef[4];
527 retval.dCoef[5] = dCoef[5] + c.dCoef[5];
528
529 return retval;
530 }
531
532 GenericAerodynamicData::GenericAerodynamicCoef
operator -(const GenericAerodynamicData::GenericAerodynamicCoef & c) const533 GenericAerodynamicData::GenericAerodynamicCoef::operator - (
534 const GenericAerodynamicData::GenericAerodynamicCoef& c) const
535 {
536 GenericAerodynamicData::GenericAerodynamicCoef retval;
537
538 retval.dCoef[0] = dCoef[0] - c.dCoef[0];
539 retval.dCoef[1] = dCoef[1] - c.dCoef[1];
540 retval.dCoef[2] = dCoef[2] - c.dCoef[2];
541 retval.dCoef[3] = dCoef[3] - c.dCoef[3];
542 retval.dCoef[4] = dCoef[4] - c.dCoef[4];
543 retval.dCoef[5] = dCoef[5] - c.dCoef[5];
544
545 return retval;
546 }
547
548 GenericAerodynamicData::GenericAerodynamicCoef
operator *(const doublereal & d) const549 GenericAerodynamicData::GenericAerodynamicCoef::operator * (const doublereal& d) const
550 {
551 GenericAerodynamicData::GenericAerodynamicCoef retval;
552
553 retval.dCoef[0] = dCoef[0]*d;
554 retval.dCoef[1] = dCoef[1]*d;
555 retval.dCoef[2] = dCoef[2]*d;
556 retval.dCoef[3] = dCoef[3]*d;
557 retval.dCoef[4] = dCoef[4]*d;
558 retval.dCoef[5] = dCoef[5]*d;
559
560 return retval;
561 }
562
563 GenericAerodynamicData::GenericAerodynamicCoef
operator /(const doublereal & d) const564 GenericAerodynamicData::GenericAerodynamicCoef::operator / (const doublereal& d) const
565 {
566 GenericAerodynamicData::GenericAerodynamicCoef retval;
567
568 retval.dCoef[0] = dCoef[0]/d;
569 retval.dCoef[1] = dCoef[1]/d;
570 retval.dCoef[2] = dCoef[2]/d;
571 retval.dCoef[3] = dCoef[3]/d;
572 retval.dCoef[4] = dCoef[4]/d;
573 retval.dCoef[5] = dCoef[5]/d;
574
575 return retval;
576 }
577
578 /* GenericAerodynamicData::GenericAerodynamicCoef - end */
579
580
581 static GenericAerodynamicData *
ReadGenericAerodynamicData(const std::string & fname,doublereal dScaleAngle,doublereal dScaleLength,bool bAlphaFirst)582 ReadGenericAerodynamicData(const std::string& fname,
583 doublereal dScaleAngle, doublereal dScaleLength, bool bAlphaFirst)
584 {
585 std::ifstream in(fname.c_str());
586
587 if (!in) {
588 silent_cerr("ReadGenericAerodynamicData: "
589 "unable to open file \"" << fname << "\""
590 << std::endl);
591 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
592 }
593
594 char buf[LINE_MAX];
595 int nAlpha, nBeta;
596 int c;
597
598 /* skip comments */
599 for (c = in.get(); c == '%' || c == '#'; c = in.get()) {
600 /* discard to end of line */
601 in.getline(buf, sizeof(buf));
602 }
603 in.putback(c);
604
605 /* get the size of the matrices */
606 in >> nAlpha >> nBeta;
607 /* discard to end of line */
608 in.getline(buf, sizeof(buf));
609
610 if (!in) {
611 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
612 "unable to read size of data matrix "
613 "from file \"" << fname << "\"" << std::endl);
614 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
615 }
616
617 if (nAlpha <= 0) {
618 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
619 "invalid number of angles of attack " << nAlpha << " "
620 "from file \"" << fname << "\"" << std::endl);
621 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
622 }
623
624 if (nBeta <= 0) {
625 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
626 "invalid number of sideslip angles " << nBeta << " "
627 "from file \"" << fname << "\"" << std::endl);
628 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
629 }
630
631 /* skip comments */
632 for (c = in.get(); c == '%' || c == '#'; c = in.get()) {
633 /* discard to end of line */
634 in.getline(buf, sizeof(buf));
635 }
636 in.putback(c);
637
638 if (!in) {
639 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
640 "unable to get to data "
641 "in file \"" << fname << "\"" << std::endl);
642 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
643 }
644
645 GenericAerodynamicData *pData = new GenericAerodynamicData;
646 pData->name = fname;
647 pData->bAlphaFirst = bAlphaFirst;
648 pData->nAlpha = nAlpha;
649 pData->nBeta = nBeta;
650
651 pData->Alpha.resize(nAlpha);
652 pData->Beta.resize(nBeta);
653 pData->Data.resize(nBeta);
654
655 doublereal dScaleForce = dScaleLength*dScaleLength;
656 doublereal dScaleMoment = dScaleForce*dScaleLength;
657 doublereal dScale[6] = {
658 dScaleForce, dScaleForce, dScaleForce,
659 dScaleMoment, dScaleMoment, dScaleMoment
660 };
661
662 /* get the matrices */
663 for (int iBeta = 0; iBeta < nBeta; iBeta++) {
664 pData->Data[iBeta].resize(nAlpha);
665
666 for (int iAlpha = 0; iAlpha < nAlpha; iAlpha++) {
667 doublereal dCoef;
668
669 /* read (and check) alpha */
670 in >> dCoef;
671 if (iBeta == 0) {
672 if (iAlpha == 0) {
673 doublereal dErr = dCoef*dScaleAngle + ::dAlphaMax[bAlphaFirst];
674 if (std::abs(dErr) > std::numeric_limits<doublereal>::epsilon()) {
675 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
676 "warning, alpha[0] != -pi/2 (error=" << 100*dErr/::dAlphaMax[bAlphaFirst] << "%)" << std::endl);
677 }
678 } else if (iAlpha == nAlpha - 1) {
679 doublereal dErr = dCoef*dScaleAngle - ::dAlphaMax[bAlphaFirst];
680 if (std::abs(dErr) > std::numeric_limits<doublereal>::epsilon()) {
681 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
682 "warning, alpha[" << iAlpha << "] != pi/2 (error=" << 100*dErr/::dAlphaMax[bAlphaFirst] << "%)" << std::endl);
683 }
684 }
685
686 pData->Alpha[iAlpha] = dCoef;
687
688 } else if (dCoef != pData->Alpha[iAlpha]) {
689 silent_cerr("ReadGenericAerodynamicData"
690 "(" << fname << "): "
691 "inconsistent data, "
692 "Alpha[" << iAlpha << "]"
693 "=" << dCoef << " "
694 "for Beta[" << iBeta << "]"
695 "=" << pData->Beta[iBeta] << " "
696 "differs from previous, "
697 << pData->Alpha[iAlpha]
698 << std::endl);
699 delete pData;
700 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
701 }
702
703 /* read (and check) beta */
704 in >> dCoef;
705 if (iAlpha == 0) {
706 if (iBeta == 0) {
707 doublereal dErr = dCoef*dScaleAngle + ::dBetaMax[bAlphaFirst];
708 if (std::abs(dErr) > std::numeric_limits<doublereal>::epsilon()) {
709 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
710 "warning, beta[0] != -pi (error=" << 100*dErr/::dBetaMax[bAlphaFirst] << "%)" << std::endl);
711 }
712 } else if (iBeta == nBeta - 1) {
713 doublereal dErr = dCoef*dScaleAngle - ::dBetaMax[bAlphaFirst];
714 if (std::abs(dErr) > std::numeric_limits<doublereal>::epsilon()) {
715 silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
716 "warning, beta[" << iBeta << "] != pi (error=" << 100*dErr/::dBetaMax[bAlphaFirst] << "%)" << std::endl);
717 }
718 }
719
720 pData->Beta[iBeta] = dCoef;
721
722 } else if (dCoef != pData->Beta[iBeta]) {
723 silent_cerr("ReadGenericAerodynamicData"
724 "(" << fname << "): "
725 "inconsistent data, "
726 "Beta[" << iBeta << "]"
727 "=" << dCoef << " "
728 "for Alpha[" << iAlpha << "]"
729 "=" << pData->Alpha[iAlpha] << " "
730 "differs from previous, "
731 << pData->Beta[iBeta]
732 << std::endl);
733 delete pData;
734 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
735 }
736
737 for (int iCoef = 0; iCoef < 6; iCoef++) {
738 in >> dCoef;
739
740 pData->Data[iBeta][iAlpha].dCoef[iCoef] = dCoef*dScale[iCoef];
741 }
742
743 /* discard to end of line */
744 if (iAlpha < nAlpha - 1 && iBeta < nBeta - 1) {
745 in.getline(buf, sizeof(buf));
746
747 if (!in) {
748 silent_cerr("ReadGenericAerodynamicData"
749 "(" << fname << "): "
750 "unable to read data past "
751 "iAlpha=" << iAlpha << ", "
752 "iBeta=" << iBeta << std::endl);
753 delete pData;
754 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
755 }
756 }
757 }
758 }
759
760 /* deg => radian */
761 for (int iAlpha = 0; iAlpha < nAlpha; iAlpha++) {
762 pData->Alpha[iAlpha] *= dScaleAngle;
763
764 if (iAlpha == 0) {
765 continue;
766 }
767
768 if ( pData->Alpha[iAlpha] <= pData->Alpha[iAlpha - 1]) {
769 silent_cerr("ReadGenericAerodynamicData"
770 "(" << fname << "): "
771 "strict ordering violated between "
772 "Alpha #" << iAlpha - 1 << " and "
773 "Alpha #" << iAlpha << std::endl);
774 delete pData;
775 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
776 }
777 }
778
779 /* deg => radian */
780 for (int iBeta = 0; iBeta < nBeta; iBeta++) {
781 pData->Beta[iBeta] *= dScaleAngle;
782
783 if (iBeta == 0) {
784 continue;
785 }
786
787 if ( pData->Beta[iBeta] <= pData->Beta[iBeta - 1]) {
788 silent_cerr("ReadGenericAerodynamicData"
789 "(" << fname << "): "
790 "strict ordering violated between "
791 "Beta #" << iBeta - 1 << " and "
792 "Beta #" << iBeta << std::endl);
793 delete pData;
794 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
795 }
796 }
797
798 return pData;
799 }
800
801 Elem *
ReadGenericAerodynamicForce(DataManager * pDM,MBDynParser & HP,const DofOwner * pDO,unsigned int uLabel)802 ReadGenericAerodynamicForce(DataManager* pDM, MBDynParser& HP,
803 const DofOwner *pDO, unsigned int uLabel)
804 {
805 /* Nodo */
806 const StructNode* pNode = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
807
808 /* The offset is in the reference frame of the node */
809 ReferenceFrame RF(pNode);
810 Vec3 f(Zero3);
811 if (HP.IsKeyWord("position")) {
812 f = HP.GetPosRel(RF);
813 }
814
815 /* the orientation is in flight mechanics (FIXME?) reference frame:
816 * X forward, Y to the right, Z down */
817 Mat3x3 Ra(Eye3);
818 if (HP.IsKeyWord("orientation")) {
819 Ra = HP.GetRotRel(RF);
820 }
821
822 /* 1. by default, which means that coefficients are only normalized
823 * by the dynamic pressure */
824 doublereal dRefSurface = 1.;
825 doublereal dRefLength = 1.;
826 bool bAlphaFirst(false);
827
828 if (HP.IsKeyWord("reference" "surface")) {
829 dRefSurface = HP.GetReal();
830 }
831
832 if (HP.IsKeyWord("reference" "length")) {
833 dRefLength = HP.GetReal();
834 }
835
836 if (HP.IsKeyWord("alpha" "first")) {
837 if (!HP.GetYesNo(bAlphaFirst)) {
838 silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
839 "invalid \"alpha first\" value at line " << HP.GetLineData() << std::endl);
840 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
841 }
842 }
843
844 /* TODO: allow to reference previously loaded data */
845 GenericAerodynamicData *pData = 0;
846 if (HP.IsKeyWord("file")) {
847 doublereal dScaleAngle = 1.;
848 if (HP.IsKeyWord("angle" "units")) {
849 if (HP.IsKeyWord("radians")) {
850 dScaleAngle = 1.;
851
852 } else if (HP.IsKeyWord("degrees")) {
853 dScaleAngle = M_PI/180.;
854
855 } else {
856 silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
857 "unknown angle unit at line " << HP.GetLineData() << std::endl);
858 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
859 }
860
861 } else if (HP.IsKeyWord("scale" "angles")) {
862 doublereal d = HP.GetReal(dScaleAngle);
863 if (d <= 0.) {
864 silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
865 "invalid angle scale factor " << d
866 << " at line " << HP.GetLineData() << std::endl);
867 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
868 }
869 dScaleAngle = d;
870 }
871
872 doublereal dScaleLength = 1.;
873 if (HP.IsKeyWord("scale" "lengths")) {
874 doublereal d = HP.GetReal(dScaleLength);
875 if (d <= 0.) {
876 silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
877 "invalid length scale factor " << d
878 << " at line " << HP.GetLineData() << std::endl);
879 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
880 }
881 dScaleLength = d;
882 }
883
884 std::string fname(HP.GetFileName());
885 pData = ReadGenericAerodynamicData(fname, dScaleAngle, dScaleLength, bAlphaFirst);
886
887 } else if (HP.IsKeyWord("reference")) {
888 silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
889 "references to generic aerodynamic data not implemented yet "
890 "at line " << HP.GetLineData() << std::endl);
891 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
892
893 } else {
894 silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
895 "keyword \"file\" or \"reference\" expected "
896 "at line " << HP.GetLineData() << std::endl);
897 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
898 }
899
900 flag fOut = pDM->fReadOutput(HP, Elem::AERODYNAMIC);
901
902 Elem *pEl = 0;
903 SAFENEWWITHCONSTRUCTOR(pEl, GenericAerodynamicForce,
904 GenericAerodynamicForce(uLabel, pDO,
905 pNode, f, Ra,
906 dRefSurface, dRefLength, bAlphaFirst,
907 pData, fOut));
908
909 return pEl;
910 }
911
912 /* GenericAerodynamicData - end */
913