1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aerodata.cc,v 1.62 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 "mynewmem.h"
35 #include "aerodata_impl.h"
36 #include "gauss.h"
37 #include "submat.h"
38 #include "aerod2.h"
39 #include "dataman.h"
40 #include "drive_.h"
41
42 /* AeroMemory - begin */
43
AeroMemory(DriveCaller * pt)44 AeroMemory::AeroMemory(DriveCaller *pt)
45 : a(0), t(0), iPoints(0), numUpdates(0), pTime(pt)
46 {
47 NO_OP;
48 }
49
~AeroMemory(void)50 AeroMemory::~AeroMemory(void)
51 {
52 if (iPoints > 0) {
53 if (pTime) {
54 SAFEDELETE(pTime);
55 }
56
57 if (a) {
58 SAFEDELETEARR(a);
59 }
60 }
61 }
62
63 #define USE_POLCOE
64 void
Predict(int i,doublereal alpha,doublereal & alf1,doublereal & alf2)65 AeroMemory::Predict(int i, doublereal alpha, doublereal &alf1, doublereal &alf2)
66 {
67 /* FIXME: this should be s, but I don't want a malloc here */
68 doublereal coe[3];
69 int s = StorageSize();
70 #ifdef USE_POLCOE
71 /*
72 * FIXME: actually this is not the order of the polynomial,
73 * but the number of coefficients, e.g. a second-order polynomial
74 * has three coefficients :)
75 */
76 integer order = s;
77 #endif /* USE_POLCOE */
78
79 ASSERT(s > 0);
80 ASSERT(i < iPoints);
81
82 doublereal *aa = a + s*i;
83 doublereal *tt = t + s*i;
84
85 aa[s-1] = alpha;
86 tt[s-1] = pTime->dGet();
87
88 if (numUpdates >= s) {
89 #ifdef USE_POLCOE
90 __FC_DECL__(polcoe)(tt, aa, &order, coe);
91 #else /* !USE_POLCOE */
92 coe[2] = ((aa[2]-aa[0])/(tt[2]-tt[0])
93 - (aa[1]-aa[0])/(tt[1]-tt[0]))/(tt[2]-tt[1]);
94 coe[1] = (aa[1]-aa[0])/(tt[1]-tt[0]) - (tt[1]+tt[0])*coe[2];
95
96 #if 0 /* not used ... */
97 coe[1] = aa[0] - tt[0]*coe[1] - tt[0]*tt[0]*coe[2];
98 #endif
99 #endif /* !USE_POLCOE */
100
101 #if 0
102 std::cerr << "aa[0:2]= " << aa[0] << "," << aa[1] << "," << aa[2] << std::endl
103 << "tt[0:2]= " << tt[0] << "," << tt[1] << "," << tt[2] << std::endl
104 << "coe[0:2]=" << coe[0] << "," << coe[1] << "," << coe[2] << std::endl;
105 #endif /* 0 */
106
107 alf1 = coe[1]+2.*coe[2]*tt[2];
108 alf2 = 2.*coe[2];
109 } else {
110 alf1 = 0.;
111 alf2 = 0.;
112 }
113 }
114
115 void
Update(int i)116 AeroMemory::Update(int i)
117 {
118 int s = StorageSize();
119 if (s > 0) {
120 /*
121 * shift back angle of attack and time
122 * for future interpolation
123 */
124 doublereal *aa = a + s*i;
125 doublereal *tt = t + s*i;
126
127 ASSERT(i < iPoints);
128
129 for (int j = 1; j < s; j++) {
130 aa[j-1] = aa[j];
131 tt[j-1] = tt[j];
132 }
133
134 if (i == 0) {
135 numUpdates++;
136 }
137 }
138 }
139
140 void
SetNumPoints(int i)141 AeroMemory::SetNumPoints(int i)
142 {
143 int s = StorageSize();
144
145 iPoints = i;
146
147 if (s > 0) {
148 SAFENEWARR(a, doublereal, 2*s*iPoints);
149 t = a + s*iPoints;
150
151 #ifdef HAVE_MEMSET
152 memset(a, 0, 2*s*iPoints*sizeof(doublereal));
153 #else /* !HAVE_MEMSET */
154 for (int i = 0; i < 2*s*iPoints; i++) {
155 a[i] = 0.;
156 }
157 #endif /* !HAVE_MEMSET */
158 }
159 }
160
161 int
GetNumPoints(void) const162 AeroMemory::GetNumPoints(void) const
163 {
164 return iPoints;
165 }
166
167 /* AeroMemory - end */
168
169 /* C81Data - begin */
170
C81Data(unsigned int uLabel)171 C81Data::C81Data(unsigned int uLabel)
172 : WithLabel(uLabel)
173 {
174 NO_OP;
175 }
176
177 /* C81Data - end */
178
179 /* AeroData - begin */
180
AeroData(int i_p,int i_dim,AeroData::UnsteadyModel u,DriveCaller * ptime)181 AeroData::AeroData(int i_p, int i_dim,
182 AeroData::UnsteadyModel u, DriveCaller *ptime)
183 : AeroMemory(ptime), unsteadyflag(u), Omega(0.)
184 {
185 // silence static analyzers
186 VAM.density = -1.;
187 VAM.sound_celerity = -1.;
188 VAM.chord = -1.;
189 VAM.force_position = 0.;
190 VAM.bc_position = 0.;
191 VAM.twist = 0.;
192
193 if (u != AeroData::STEADY) {
194 if (ptime == 0) {
195 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
196 }
197 }
198
199 SetNumPoints(i_p*i_dim);
200 }
201
~AeroData(void)202 AeroData::~AeroData(void)
203 {
204 NO_OP;
205 }
206
207 AeroData::UnsteadyModel
Unsteady(void) const208 AeroData::Unsteady(void) const
209 {
210 return unsteadyflag;
211 }
212
213 void
SetAirData(const doublereal & rho,const doublereal & c)214 AeroData::SetAirData(const doublereal& rho, const doublereal& c)
215 {
216 VAM.density = rho;
217 VAM.sound_celerity = c;
218 }
219
220 int
StorageSize(void) const221 AeroData::StorageSize(void) const
222 {
223 switch (unsteadyflag) {
224 case AeroData::HARRIS:
225 case AeroData::BIELAWA:
226 return 3;
227
228 case AeroData::STEADY:
229 return 0;
230
231 default:
232 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
233 }
234 }
235
236 void
SetSectionData(const doublereal & abscissa,const doublereal & chord,const doublereal & forcepoint,const doublereal & velocitypoint,const doublereal & twist,const doublereal & omega)237 AeroData::SetSectionData(const doublereal& abscissa,
238 const doublereal& chord,
239 const doublereal& forcepoint,
240 const doublereal& velocitypoint,
241 const doublereal& twist,
242 const doublereal& omega)
243 {
244 VAM.chord = chord;
245 VAM.force_position = forcepoint;
246 VAM.bc_position = velocitypoint;
247 VAM.twist = twist;
248 Omega = omega;
249 }
250
251 std::ostream&
RestartUnsteady(std::ostream & out) const252 AeroData::RestartUnsteady(std::ostream& out) const
253 {
254 switch (unsteadyflag) {
255 case AeroData::HARRIS:
256 out << ", unsteady, harris";
257 break;
258
259 case AeroData::BIELAWA:
260 out << ", unsteady, bielawa";
261 break;
262
263 default:
264 break;
265 }
266
267 return out;
268 }
269
270 int
GetForcesJacForwardDiff_int(int i,const doublereal * W,doublereal * TNG0,Mat6x6 & J,outa_t & OUTA)271 AeroData::GetForcesJacForwardDiff_int(int i, const doublereal* W, doublereal* TNG0, Mat6x6& J, outa_t& OUTA)
272 {
273 const doublereal epsilon = 1.e-3;
274 const doublereal nu = 1.e-9;
275
276 Vec3 V(&W[0]);
277 Vec3 Omega(&W[3]);
278
279 doublereal dv = V.Norm();
280 doublereal dw = Omega.Norm();
281
282 doublereal TNG[6];
283
284 GetForces(i, W, TNG0, OUTA);
285
286 doublereal *WW = const_cast<doublereal *>(W);
287 for (unsigned int iColm1 = 0; iColm1 < 6; iColm1++) {
288 doublereal dorig;
289 doublereal delta;
290
291 dorig = WW[iColm1];
292 if (iColm1 < 3) {
293 delta = dv*epsilon + nu;
294
295 } else {
296 delta = dw*epsilon + nu;
297 }
298 WW[iColm1] = dorig + delta;
299
300 GetForces(i, WW, TNG, OUTA);
301
302 for (unsigned int iRowm1 = 0; iRowm1 < 6; iRowm1++) {
303 J.Put(iRowm1 + 1, iColm1 + 1, (TNG[iRowm1] - TNG0[iRowm1])/delta);
304 }
305
306 WW[iColm1] = dorig;
307 }
308
309 return 0;
310 }
311
312 int
GetForcesJacCenteredDiff_int(int i,const doublereal * W,doublereal * TNG0,Mat6x6 & J,outa_t & OUTA)313 AeroData::GetForcesJacCenteredDiff_int(int i, const doublereal* W, doublereal* TNG0, Mat6x6& J, outa_t& OUTA)
314 {
315 const doublereal epsilon = 1.e-3;
316 const doublereal nu = 1.e-9;
317
318 Vec3 V(&W[0]);
319 Vec3 Omega(&W[3]);
320
321 doublereal dv = V.Norm();
322 doublereal dw = Omega.Norm();
323
324 doublereal TNGp[6];
325 doublereal TNGm[6];
326
327 GetForces(i, W, TNG0, OUTA);
328
329 doublereal *WW = const_cast<doublereal *>(W);
330 for (unsigned int iColm1 = 0; iColm1 < 6; iColm1++) {
331 doublereal dorig;
332 doublereal delta;
333
334 dorig = WW[iColm1];
335 if (iColm1 < 3) {
336 delta = dv*epsilon + nu;
337
338 } else {
339 delta = dw*epsilon + nu;
340 }
341
342 WW[iColm1] = dorig + delta;
343
344 GetForces(i, WW, TNGp, OUTA);
345
346 WW[iColm1] = dorig - delta;
347
348 GetForces(i, WW, TNGm, OUTA);
349
350 doublereal delta2 = 2.*delta;
351 for (unsigned int iRowm1 = 0; iRowm1 < 6; iRowm1++) {
352 J.Put(iRowm1 + 1, iColm1 + 1, (TNGp[iRowm1] - TNGm[iRowm1])/delta2);
353 }
354
355 WW[iColm1] = dorig;
356 }
357
358 return 0;
359 }
360
361 unsigned int
iGetNumDof(void) const362 AeroData::iGetNumDof(void) const
363 {
364 return 0;
365 }
366
367 DofOrder::Order
GetDofType(unsigned int) const368 AeroData::GetDofType(unsigned int) const
369 {
370 silent_cerr("AeroData: "
371 "GetDofType() is undefined because aerodynamic data "
372 "has no degrees of freedom" << std::endl);
373 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
374 }
375
376 int
GetForces(int i,const doublereal * W,doublereal * TNG,outa_t & OUTA)377 AeroData::GetForces(int i, const doublereal* W, doublereal* TNG, outa_t& OUTA)
378 {
379 silent_cerr("AeroData: GetForces() is undefined" << std::endl);
380 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
381 }
382
383 int
GetForcesJac(int i,const doublereal * W,doublereal * TNG,Mat6x6 & J,outa_t & OUTA)384 AeroData::GetForcesJac(int i, const doublereal* W, doublereal* TNG, Mat6x6& J, outa_t& OUTA)
385 {
386 silent_cerr("AeroData: GetForcesJac() is undefined" << std::endl);
387 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
388 }
389
390 void
AssRes(SubVectorHandler & WorkVec,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr,integer iFirstIndex,integer iFirstSubIndex,int i,const doublereal * W,doublereal * TNG,outa_t & OUTA)391 AeroData::AssRes(SubVectorHandler& WorkVec,
392 doublereal dCoef,
393 const VectorHandler& XCurr,
394 const VectorHandler& XPrimeCurr,
395 integer iFirstIndex, integer iFirstSubIndex,
396 int i, const doublereal* W, doublereal* TNG, outa_t& OUTA)
397 {
398 silent_cerr("AeroData: AssRes() is undefined" << std::endl);
399 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
400 }
401
402 void
AssJac(FullSubMatrixHandler & WorkMat,doublereal dCoef,const VectorHandler & XCurr,const VectorHandler & XPrimeCurr,integer iFirstIndex,integer iFirstSubIndex,const Mat3xN & vx,const Mat3xN & wx,Mat3xN & fq,Mat3xN & cq,int i,const doublereal * W,doublereal * TNG,Mat6x6 & J,outa_t & OUTA)403 AeroData::AssJac(FullSubMatrixHandler& WorkMat,
404 doublereal dCoef,
405 const VectorHandler& XCurr,
406 const VectorHandler& XPrimeCurr,
407 integer iFirstIndex, integer iFirstSubIndex,
408 const Mat3xN& vx, const Mat3xN& wx, Mat3xN& fq, Mat3xN& cq,
409 int i, const doublereal* W, doublereal* TNG, Mat6x6& J, outa_t& OUTA)
410 {
411 silent_cerr("AeroData: AssJac() is undefined" << std::endl);
412 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
413 }
414
415 void
AfterConvergence(int i,const VectorHandler & X,const VectorHandler & XP)416 AeroData::AfterConvergence(
417 int i,
418 const VectorHandler& X,
419 const VectorHandler& XP)
420 {
421 NO_OP;
422 }
423
424 /* AeroData - end */
425
426 static AeroData::UnsteadyModel
ReadUnsteadyFlag(MBDynParser & HP)427 ReadUnsteadyFlag(MBDynParser& HP)
428 {
429 if (HP.IsArg()) {
430 AeroData::UnsteadyModel eInst = AeroData::STEADY;
431 if (HP.IsKeyWord("unsteady")) {
432 /*
433 * swallow "unsteady" keyword
434 */
435 if (HP.IsKeyWord("steady")) {
436 eInst = AeroData::STEADY;
437 } else if (HP.IsKeyWord("harris")) {
438 eInst = AeroData::HARRIS;
439 } else if (HP.IsKeyWord("bielawa")) {
440 eInst = AeroData::BIELAWA;
441 } else {
442 /* demote to pedantic, because the integer
443 * form allows to change unsteady model
444 * parametrically (while waiting for string
445 * vars) */
446 pedantic_cerr("deprecated unsteady model "
447 "given by integer number;"
448 " use \"steady\", \"Harris\" or \"Bielawa\" "
449 "instead, at line " << HP.GetLineData()
450 << std::endl);
451
452 int i = HP.GetInt();
453 if (i < AeroData::STEADY || i >= AeroData::LAST) {
454 silent_cerr("illegal unsteady flag "
455 "numeric value " << i
456 << " at line "
457 << HP.GetLineData()
458 << std::endl);
459 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
460 }
461 eInst = AeroData::UnsteadyModel(i);
462 }
463 }
464
465 switch (eInst) {
466 case AeroData::STEADY:
467 case AeroData::BIELAWA:
468 break;
469
470 case AeroData::HARRIS:
471 silent_cerr("\"Harris\" unsteady aerodynamics "
472 "are not available at line "
473 << HP.GetLineData()
474 << std::endl);
475 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
476
477 default:
478 silent_cerr("illegal unsteady flag at line "
479 << HP.GetLineData() << std::endl);
480 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
481 }
482
483 /*
484 * unsteady flag
485 */
486 return eInst;
487 }
488
489 /*
490 * default: no unsteady ...
491 */
492 return AeroData::STEADY;
493 }
494
495 static void
ReadC81MultipleAeroData(DataManager * pDM,MBDynParser & HP,AeroData ** aerodata,int iGP,int iDim,bool bInterp=false)496 ReadC81MultipleAeroData(DataManager* pDM, MBDynParser& HP, AeroData** aerodata,
497 int iGP, int iDim, bool bInterp = false)
498 {
499 doublereal dcltol = 1.e-6;
500 if (bInterp && HP.IsKeyWord("tolerance")) {
501 dcltol = HP.GetReal();
502 if (dcltol <= 0.) {
503 silent_cerr("ReadC81MultipleAeroData: "
504 "invalid c81 data tolerance at line "
505 << HP.GetLineData() << std::endl);
506 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
507 }
508 }
509
510 integer nProfiles = HP.GetInt();
511 if (nProfiles <= 0) {
512 silent_cerr("Need at least one airfoil at line "
513 << HP.GetLineData() << std::endl);
514 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
515 }
516
517 std::vector<unsigned> profiles(nProfiles);
518 std::vector<doublereal> upper_bounds(nProfiles);
519 std::vector<const c81_data *> data(nProfiles);
520
521 for (int i = 0; i < nProfiles; i++) {
522 profiles[i] = (unsigned)HP.GetInt();
523 upper_bounds[i] = HP.GetReal();
524 if (upper_bounds[i] <= -1.) {
525 if (!bInterp || upper_bounds[i] < -1.) {
526 silent_cerr("upper bound "
527 << i + 1 << " = " << upper_bounds[i]
528 << " too small at line "
529 << HP.GetLineData()
530 << std::endl);
531 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
532 }
533
534 } else if (upper_bounds[i] > 1.) {
535 silent_cerr("upper bound "
536 << i + 1 << " = "
537 << upper_bounds[i]
538 << " too large at line "
539 << HP.GetLineData()
540 << std::endl);
541 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
542
543 } else if (i > 0 && upper_bounds[i] <= upper_bounds[i-1]) {
544 silent_cerr("upper bound "
545 << i+1 << " = " << upper_bounds[i]
546 << " not in increasing order "
547 "at line " << HP.GetLineData()
548 << std::endl);
549 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
550 }
551 data[i] = HP.GetC81Data(profiles[i]);
552 if (data[i] == 0) {
553 silent_cerr("Unable to find airfoil "
554 << profiles[i] << " at line "
555 << HP.GetLineData() << std::endl);
556 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
557 }
558 DEBUGLCOUT(MYDEBUG_INPUT, "airfoil data " << i+1
559 << " is from file c81 " << profiles[i]
560 << std::endl);
561 }
562
563 if (upper_bounds[nProfiles - 1] != 1.) {
564 silent_cerr("warning: the last upper bound should be 1.0 "
565 "at line " << HP.GetLineData() << std::endl);
566 }
567
568 AeroData::UnsteadyModel eInst = ReadUnsteadyFlag(HP);
569 DriveCaller *ptime = 0;
570 if (eInst != AeroData::STEADY) {
571 SAFENEWWITHCONSTRUCTOR(ptime,
572 TimeDriveCaller,
573 TimeDriveCaller(pDM->pGetDrvHdl()));
574 }
575
576 if (bInterp) {
577 SAFENEWWITHCONSTRUCTOR(*aerodata,
578 C81InterpolatedAeroData,
579 C81InterpolatedAeroData(iGP, iDim,
580 eInst, profiles, upper_bounds, data,
581 dcltol, ptime));
582
583 } else {
584 SAFENEWWITHCONSTRUCTOR(*aerodata,
585 C81MultipleAeroData,
586 C81MultipleAeroData(iGP, iDim,
587 eInst, profiles, upper_bounds, data,
588 ptime));
589 }
590 }
591
592 void
ReadAeroData(DataManager * pDM,MBDynParser & HP,int iDim,Shape ** ppChord,Shape ** ppForce,Shape ** ppVelocity,Shape ** ppTwist,Shape ** ppTipLoss,integer * piNumber,DriveCaller ** ppDC,AeroData ** aerodata)593 ReadAeroData(DataManager* pDM, MBDynParser& HP, int iDim,
594 Shape** ppChord, Shape** ppForce,
595 Shape** ppVelocity, Shape** ppTwist, Shape** ppTipLoss,
596 integer* piNumber, DriveCaller** ppDC,
597 AeroData** aerodata)
598 {
599 DEBUGCOUTFNAME("ReadAeroData");
600
601 /* Keywords */
602 const char* sKeyWords[] = {
603 "naca0012",
604 "rae9671",
605 "c81",
606
607 "theodorsen",
608
609 0
610 };
611
612 /* enum delle parole chiave */
613 enum KeyWords {
614 UNKNOWN = -1,
615 NACA0012 = 0,
616 RAE9671,
617 C81,
618
619 THEODORSEN,
620
621 LASTKEYWORD
622 };
623
624 /* tabella delle parole chiave */
625 KeyTable K(HP, sKeyWords);
626
627 *ppChord = ReadShape(HP);
628 *ppForce = ReadShape(HP);
629 *ppVelocity = ReadShape(HP);
630 *ppTwist = ReadShape(HP);
631
632 *ppTipLoss = 0;
633 if (HP.IsKeyWord("tip" "loss")) {
634 *ppTipLoss = ReadShape(HP);
635
636 } else {
637 SAFENEWWITHCONSTRUCTOR(*ppTipLoss, ConstShape1D,
638 ConstShape1D(1.));
639 }
640
641 *piNumber = HP.GetInt();
642 if (*piNumber <= 0) {
643 silent_cerr("need at least 1 Gauss integration point at line "
644 << HP.GetLineData() << std::endl);
645 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
646 }
647 DEBUGLCOUT(MYDEBUG_INPUT, "Gauss points number: "
648 << *piNumber << std::endl);
649
650 if (HP.IsKeyWord("control")) {
651 /* Driver di un'eventuale controllo */
652 *ppDC = HP.GetDriveCaller();
653
654 } else {
655 SAFENEW(*ppDC, NullDriveCaller);
656 }
657
658 if (HP.IsArg()) {
659 switch (HP.IsKeyWord()) {
660 default:
661 silent_cerr("unknown airfoil type \"" << HP.GetString()
662 << "\" at line " << HP.GetLineData()
663 << "; using default (NACA0012)"
664 << std::endl);
665
666 case NACA0012: {
667 DEBUGLCOUT(MYDEBUG_INPUT,
668 "airfoil is NACA0012" << std::endl);
669
670 AeroData::UnsteadyModel eInst = ReadUnsteadyFlag(HP);
671 DriveCaller *ptime = 0;
672 if (eInst != AeroData::STEADY) {
673 SAFENEWWITHCONSTRUCTOR(ptime, TimeDriveCaller,
674 TimeDriveCaller(pDM->pGetDrvHdl()));
675 }
676 SAFENEWWITHCONSTRUCTOR(*aerodata,
677 STAHRAeroData,
678 STAHRAeroData(*piNumber, iDim, eInst, 1, ptime));
679 break;
680 }
681
682 case RAE9671: {
683 DEBUGLCOUT(MYDEBUG_INPUT,
684 "airfoil is RAE9671" << std::endl);
685
686 AeroData::UnsteadyModel eInst = ReadUnsteadyFlag(HP);
687 DriveCaller *ptime = 0;
688 if (eInst != AeroData::STEADY) {
689 SAFENEWWITHCONSTRUCTOR(ptime, TimeDriveCaller,
690 TimeDriveCaller(pDM->pGetDrvHdl()));
691 }
692 SAFENEWWITHCONSTRUCTOR(*aerodata,
693 STAHRAeroData,
694 STAHRAeroData(*piNumber, iDim, eInst, 2, ptime));
695 break;
696 }
697
698 /*
699 * uso tabelle standard, che vengono cercate in base all'indice
700 * (sistemare)
701 */
702 case C81:
703 if (HP.IsKeyWord("multiple")) {
704 ReadC81MultipleAeroData(pDM, HP, aerodata, *piNumber, iDim);
705
706 } else if (HP.IsKeyWord("interpolated")) {
707 ReadC81MultipleAeroData(pDM, HP, aerodata, *piNumber, iDim, true);
708
709 } else {
710 unsigned iProfile = (unsigned)HP.GetInt();
711 const c81_data* data = HP.GetC81Data(iProfile);
712
713 DEBUGLCOUT(MYDEBUG_INPUT,
714 "airfoil data is from file c81 "
715 << iProfile << std::endl);
716 AeroData::UnsteadyModel
717 eInst = ReadUnsteadyFlag(HP);
718 DriveCaller *ptime = 0;
719 if (eInst != AeroData::STEADY) {
720 SAFENEWWITHCONSTRUCTOR(ptime,
721 TimeDriveCaller,
722 TimeDriveCaller(pDM->pGetDrvHdl()));
723 }
724 SAFENEWWITHCONSTRUCTOR(*aerodata,
725 C81AeroData,
726 C81AeroData(*piNumber, iDim,
727 eInst, iProfile,
728 data, ptime));
729 }
730 break;
731
732 case THEODORSEN: {
733 if (!HP.IsKeyWord("c81")) {
734 silent_cerr("Theodorsen aerodata: warning, assuming \"c81\" at line "
735 << HP.GetLineData() << std::endl);
736 }
737
738 AeroData *pa = 0;
739 if (HP.IsKeyWord("multiple")) {
740 ReadC81MultipleAeroData(pDM, HP, &pa, *piNumber, iDim);
741
742 } else if (HP.IsKeyWord("interpolated")) {
743 ReadC81MultipleAeroData(pDM, HP, &pa, *piNumber, iDim, true);
744
745 } else {
746 unsigned iProfile = (unsigned)HP.GetInt();
747 const c81_data* data = HP.GetC81Data(iProfile);
748
749 DEBUGLCOUT(MYDEBUG_INPUT,
750 "airfoil data is from file c81 "
751 << iProfile << std::endl);
752 AeroData::UnsteadyModel
753 eInst = ReadUnsteadyFlag(HP);
754 DriveCaller *ptime = 0;
755 if (eInst != AeroData::STEADY) {
756 SAFENEWWITHCONSTRUCTOR(ptime,
757 TimeDriveCaller,
758 TimeDriveCaller(pDM->pGetDrvHdl()));
759 }
760 SAFENEWWITHCONSTRUCTOR(pa,
761 C81AeroData,
762 C81AeroData(*piNumber, iDim,
763 eInst, iProfile,
764 data, ptime));
765 }
766
767 DriveCaller *ptime = 0;
768 SAFENEWWITHCONSTRUCTOR(ptime,
769 TimeDriveCaller,
770 TimeDriveCaller(pDM->pGetDrvHdl()));
771
772 SAFENEWWITHCONSTRUCTOR(*aerodata,
773 TheodorsenAeroData,
774 TheodorsenAeroData(*piNumber, iDim,
775 pa, ptime));
776 } break;
777 }
778
779 } else {
780 /* FIXME: better abort! */
781 silent_cerr("missing airfoil type at line "
782 << HP.GetLineData()
783 << "; using default (NACA0012)"
784 << std::endl);
785
786 AeroData::UnsteadyModel eInst = ReadUnsteadyFlag(HP);
787 SAFENEWWITHCONSTRUCTOR(*aerodata,
788 STAHRAeroData, STAHRAeroData(*piNumber, iDim, eInst, 1));
789 }
790 }
791
792