1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aerodyn.cc,v 1.51 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
33  *
34  * - Proprieta' dell'aria:
35  * elemento unico, contiene:
36  *   - direzione e modulo del vento relativo. Il modulo e' associato
37  *     ad un driver che ne consente la variazione, per consentire la
38  *     simulazione di un transitorio di galleria.
39  *     Verra' aggiunta una variazione della direzione per simulare la raffica.
40  *   - densita' dell'aria
41  *   - celerita' del suono
42  *   - ...
43  *
44  * - Elemento di rotore:
45  * classe virtuale che possiede alcuni dati topologici e geometrici di rotore
46  * ed i metodi per calcolare grandezze utili agli elementi aerodinamici che
47  * fanno parte di un rotore
48  *
49  * - Elemento aerodinamico:
50  * stazionario o quasi-stazionario, con metodo p-k di ordine 0, 1 o 2,
51  * associato ad un corpo rigido, basato sulla strip theory.
52  *
53  * - Elemento aerodinamico:
54  * analogo al precedente, ma associato alla trave a Volumi Finiti a tre nodi
55  *
56  * - Elemento aerodinamico instazionario:
57  * in fase di sviluppo, modella dinamicamente alcuni stati a dare
58  * il comportamento instazionario di una superficie aerodinamica
59  * modellata con la strip theory
60  *
61  */
62 
63 #include "mbconfig.h"           /* This goes first in every *.c,*.cc file */
64 
65 #include "aerodyn_.h"
66 #include "dataman.h"
67 #include "tpldrive.h"
68 
69 #include "gust.h"
70 
71 /* AirProperties - begin */
72 
AirProperties(const TplDriveCaller<Vec3> * pDC,std::vector<const Gust * > & g,const RigidBodyKinematics * pRBK,flag fOut)73 AirProperties::AirProperties(const TplDriveCaller<Vec3>* pDC,
74 		std::vector<const Gust *>& g,
75 		const RigidBodyKinematics *pRBK,
76 		flag fOut)
77 : Elem(1, fOut),
78 InitialAssemblyElem(1, fOut),
79 TplDriveOwner<Vec3>(pDC),
80 Velocity(Zero3),
81 gust(g),
82 pRBK(pRBK)
83 {
84 	for (std::vector<const Gust *>::iterator i = gust.begin();
85 		i != gust.end(); ++i)
86 	{
87 		const_cast<Gust *>(*i)->SetAirProperties(this);
88 	}
89 }
90 
~AirProperties(void)91 AirProperties::~AirProperties(void)
92 {
93 	for (std::vector<const Gust *>::iterator i = gust.begin();
94 		i != gust.end(); ++i)
95 	{
96 		SAFEDELETE(*i);
97 	}
98 }
99 
100 void
AddGust(const Gust * pG)101 AirProperties::AddGust(const Gust *pG)
102 {
103 	gust.insert(gust.end(), pG);
104 	const_cast<Gust *>(pG)->SetAirProperties(this);
105 }
106 
107 /* Scrive il contributo dell'elemento al file di restart */
108 std::ostream&
Restart(std::ostream & out) const109 AirProperties::Restart(std::ostream& out) const
110 {
111 	pGetDriveCaller()->Restart(out);
112 	if (!gust.empty()) {
113 		for (std::vector<const Gust *>::const_iterator i = gust.begin();
114 			i != gust.end(); ++i)
115 		{
116 			out << ", ", (*i)->Restart(out);
117 		}
118 	}
119 	return out << ';' << std::endl;
120 }
121 
122 /* Dimensioni del workspace */
123 void
WorkSpaceDim(integer * piNumRows,integer * piNumCols) const124 AirProperties::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
125 {
126 	*piNumRows = 0;
127 	*piNumCols = 0;
128 }
129 
130 /* assemblaggio jacobiano */
131 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal,const VectorHandler &,const VectorHandler &)132 AirProperties::AssJac(VariableSubMatrixHandler& WorkMat,
133 		doublereal /* dCoef */ ,
134 		const VectorHandler& /* XCurr */ ,
135 		const VectorHandler& /* XPrimeCurr */ )
136 {
137 	WorkMat.SetNullMatrix();
138 	return WorkMat;
139 }
140 
141 /* assemblaggio residuo */
142 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal,const VectorHandler &,const VectorHandler &)143 AirProperties::AssRes(SubVectorHandler& WorkVec,
144 		doublereal /* dCoef */ ,
145 		const VectorHandler& /* XCurr */ ,
146 		const VectorHandler& /* XPrimeCurr */ )
147 {
148   	WorkVec.Resize(0);
149 
150 	/* Approfitto del fatto che AirProperties viene aggiornato prima
151 	 * degli altri elementi (vedi l'enum Elem::Type e la sequenza di
152 	 * assemblaggio) per fargli calcolare Velocity una volta per tutte.
153 	 * Quindi, quando viene chiamata GetVelocity(void),
154 	 * questa restituisce un reference alla velocita' con il
155 	 * minimo overhead */
156 	Velocity = Get();
157 
158 	return WorkVec;
159 }
160 
161 /* Numero di GDL iniziali */
162 unsigned int
iGetInitialNumDof(void) const163 AirProperties::iGetInitialNumDof(void) const
164 {
165 	return 0;
166 }
167 
168 /* Dimensioni initiali del workspace */
169 void
InitialWorkSpaceDim(integer * piNumRows,integer * piNumCols) const170 AirProperties::InitialWorkSpaceDim(integer* piNumRows,
171 		integer* piNumCols) const
172 {
173 	*piNumRows = 0;
174 	*piNumCols = 0;
175 }
176 
177 /* assemblaggio jacobiano */
178 VariableSubMatrixHandler&
InitialAssJac(VariableSubMatrixHandler & WorkMat,const VectorHandler &)179 AirProperties::InitialAssJac(VariableSubMatrixHandler& WorkMat,
180 		const VectorHandler& /* XCurr */ )
181 {
182 	WorkMat.SetNullMatrix();
183 	return WorkMat;
184 }
185 
186 /* assemblaggio residuo */
187 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)188 AirProperties::InitialAssRes(SubVectorHandler& WorkVec,
189 		const VectorHandler& XCurr)
190 {
191 	/* Chiama AssRes, che si limita ad aggiornare la velocita' */
192 	return AssRes(WorkVec, 1, XCurr, XCurr);
193 }
194 
195 void
Output(OutputHandler & OH) const196 AirProperties::Output(OutputHandler& OH) const
197 {
198 	if (bToBeOutput()) {
199 		OH.AirProps() << std::setw(8) << GetLabel()
200 			<< " " << dGetAirDensity(Zero3)
201 			<< " " << dGetSoundSpeed(Zero3)
202 			<< " " << GetVelocity(Zero3)
203 			<< std::endl;
204 	}
205 }
206 
207 Vec3
GetVelocity(const Vec3 & X) const208 AirProperties::GetVelocity(const Vec3& X) const
209 {
210 	Vec3 V;
211 
212 	GetVelocity(X, V);
213 
214 	return V;
215 }
216 
217 bool
GetVelocity(const Vec3 & X,Vec3 & V) const218 AirProperties::GetVelocity(const Vec3& X, Vec3& V) const
219 {
220 	Vec3 Xabs;
221 	if (pRBK) {
222 		// X is the position of the point in the relative frame
223 		// Xabs is the position of the point in the absolute frame
224 		Vec3 Xabs = pRBK->GetX();
225 		Xabs += pRBK->GetR()*X;
226 
227 	} else {
228 		Xabs = X;
229 	}
230 
231 	V = Velocity;
232 	for (std::vector<const Gust *>::const_iterator i = gust.begin();
233 		i != gust.end(); ++i)
234 	{
235 		Vec3 VV;
236 		if ((*i)->GetVelocity(Xabs, VV)) {
237 			V += VV;
238 		}
239 	}
240 
241 	if (pRBK) {
242 		// V is the velocity of the point in the absolute frame,
243 		// projected in the relative reference frame,
244 		// plus the airstream speed, minus the relative frame velocity
245 
246 		// the velocity is projected in the relative frame
247 		V = pRBK->GetR().MulTV(V);
248 
249 		// the reference velocity is subtracted
250 		V -= pRBK->GetV();
251 
252 		// the contribution of the angular velocity is subtracted
253 		V -= pRBK->GetW().Cross(X);
254 	}
255 
256 	return true;
257 }
258 
259 /* Dati privati */
260 unsigned int
iGetNumPrivData(void) const261 AirProperties::iGetNumPrivData(void) const
262 {
263 	/* 3 components + module */
264 	return 4;
265 }
266 
267 unsigned int
iGetPrivDataIdx(const char * s) const268 AirProperties::iGetPrivDataIdx(const char *s) const
269 {
270 	ASSERT(s != NULL);
271 
272 	if (strcasecmp(s, "vxinf") == 0) {
273 		return 1;
274 
275 	} else if (strcasecmp(s, "vyinf") == 0) {
276 		return 2;
277 
278 	} else if (strcasecmp(s, "vzinf") == 0) {
279 		return 3;
280 
281 	} else if (strcasecmp(s, "vinf") == 0) {
282 		return 4;
283 	}
284 
285 	return 0;
286 }
287 
288 doublereal
dGetPrivData(unsigned int i) const289 AirProperties::dGetPrivData(unsigned int i) const
290 {
291 	switch (i) {
292 	case 1:
293 	case 2:
294 	case 3:
295 		return Velocity(i);
296 
297 	case 4:
298 		return Velocity.Norm();
299 
300 	default:
301 		silent_cerr("AirProperties(" << GetLabel() << "): "
302 			"illegal property " << i << std::endl);
303 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
304 	}
305 }
306 
307 /* AirProperties - end */
308 
309 
310 /* BasicAirProperties - begin */
311 
BasicAirProperties(const TplDriveCaller<Vec3> * pDC,const DriveCaller * pRho,doublereal dSS,std::vector<const Gust * > & g,const RigidBodyKinematics * pRBK,flag fOut)312 BasicAirProperties::BasicAirProperties(const TplDriveCaller<Vec3>* pDC,
313 	const DriveCaller *pRho, doublereal dSS, std::vector<const Gust *>& g,
314 	const RigidBodyKinematics *pRBK, flag fOut)
315 : Elem(1, fOut),
316 AirProperties(pDC, g, pRBK, fOut),
317 AirDensity(pRho),
318 dSoundSpeed(dSS)
319 {
320 	NO_OP;
321 }
322 
~BasicAirProperties(void)323 BasicAirProperties::~BasicAirProperties(void)
324 {
325 	NO_OP;
326 }
327 
328 /* Scrive il contributo dell'elemento al file di restart */
329 std::ostream&
Restart(std::ostream & out) const330 BasicAirProperties::Restart(std::ostream& out) const
331 {
332 	out << "  air properties: ",
333 		AirDensity.pGetDriveCaller()->Restart(out) << ", " << dSoundSpeed << ", ";
334 	return AirProperties::Restart(out);
335 }
336 
337 doublereal
dGetAirDensity(const Vec3 &) const338 BasicAirProperties::dGetAirDensity(const Vec3& /* X */ ) const
339 {
340 	return AirDensity.dGet();
341 }
342 
343 doublereal
dGetAirPressure(const Vec3 &) const344 BasicAirProperties::dGetAirPressure(const Vec3& /* X */ ) const
345 {
346 	return -1.;
347 }
348 
349 doublereal
dGetAirTemperature(const Vec3 &) const350 BasicAirProperties::dGetAirTemperature(const Vec3& /* X */ ) const
351 {
352 	return -1.;
353 }
354 
355 doublereal
dGetSoundSpeed(const Vec3 &) const356 BasicAirProperties::dGetSoundSpeed(const Vec3& /* X */ ) const
357 {
358 	return dSoundSpeed;
359 }
360 
361 bool
GetAirProps(const Vec3 & X,doublereal & rho,doublereal & c,doublereal & p,doublereal & T) const362 BasicAirProperties::GetAirProps(const Vec3& X, doublereal& rho,
363 		doublereal& c, doublereal& p, doublereal& T) const
364 {
365 	/* FIXME */
366 	rho = dGetAirDensity(X);
367 	c = dGetSoundSpeed(X);
368 	p = -1.;
369 	T = -1.;
370 
371 	return true;
372 }
373 
374 /* BasicAirProperties - end */
375 
376 
377 /* StdAirProperties - begin */
378 
StdAirProperties(const TplDriveCaller<Vec3> * pDC,doublereal PRef,const DriveCaller * RhoRef,doublereal TRef,doublereal a,doublereal R,doublereal g0,doublereal z0,doublereal z1,doublereal z2,std::vector<const Gust * > & g,const RigidBodyKinematics * pRBK,flag fOut)379 StdAirProperties::StdAirProperties(const TplDriveCaller<Vec3>* pDC,
380 	doublereal PRef, const DriveCaller *RhoRef, doublereal TRef,
381 	doublereal a, doublereal R, doublereal g0,
382 	doublereal z0, doublereal z1, doublereal z2,
383 	std::vector<const Gust *>& g, const RigidBodyKinematics *pRBK, flag fOut)
384 : Elem(1, fOut),
385 AirProperties(pDC, g, pRBK, fOut),
386 PRef(PRef),
387 RhoRef(RhoRef),
388 TRef(TRef),
389 a(a),
390 R(R),
391 g0(g0),
392 z0(z0),
393 z1(z1),
394 z2(z2)
395 {
396 	ASSERT(PRef > 0.);
397 	ASSERT(RhoRef != NULL);
398 	ASSERT(TRef > 0.);
399 	ASSERT(R > 0.);
400 	ASSERT(g0 > 0.);
401 	ASSERT(z1 > 0.);
402 	ASSERT(z2 > z1);
403 }
404 
~StdAirProperties(void)405 StdAirProperties::~StdAirProperties(void)
406 {
407 	delete RhoRef;
408 }
409 
410 /* Scrive il contributo dell'elemento al file di restart */
411 std::ostream&
Restart(std::ostream & out) const412 StdAirProperties::Restart(std::ostream& out) const
413 {
414 	out << "  air properties: std, "
415 		<< PRef << ", ",
416 		RhoRef->Restart(out) << ", "
417 		<< TRef << ", "
418 		<< a << ", "
419 		<< R << ", "
420 		<< g0 << ", "
421 		<< z1 << ", "
422 		<< z2 << ", ";
423 	if (z0 != 0.) {
424 		out << "reference altitude, " << z0 << ", ";
425 	}
426 	return AirProperties::Restart(out);
427 }
428 
429 doublereal
dGetAirDensity(const Vec3 & X) const430 StdAirProperties::dGetAirDensity(const Vec3& X) const
431 {
432 	doublereal rho, c, p, T;
433 
434 	GetAirProps(X, rho, c, p, T);
435 
436 	return rho;
437 }
438 
439 doublereal
dGetAirPressure(const Vec3 & X) const440 StdAirProperties::dGetAirPressure(const Vec3& X) const
441 {
442 	doublereal rho, c, p, T;
443 
444 	GetAirProps(X, rho, c, p, T);
445 
446 	return p;
447 }
448 
449 doublereal
dGetAirTemperature(const Vec3 & X) const450 StdAirProperties::dGetAirTemperature(const Vec3& X) const
451 {
452 	doublereal rho, c, p, T;
453 
454 	GetAirProps(X, rho, c, p, T);
455 
456 	return T;
457 }
458 
459 doublereal
dGetSoundSpeed(const Vec3 & X) const460 StdAirProperties::dGetSoundSpeed(const Vec3& X) const
461 {
462 	doublereal rho, c, p, T;
463 
464 	GetAirProps(X, rho, c, p, T);
465 
466 	return c;
467 }
468 
469 bool
GetAirProps(const Vec3 & X,doublereal & rho,doublereal & c,doublereal & p,doublereal & T) const470 StdAirProperties::GetAirProps(const Vec3& X, doublereal& rho,
471 		doublereal& c, doublereal& p, doublereal& T) const
472 {
473 	/* FIXME */
474 	doublereal z = X.dGet(3) + z0;
475 
476 	doublereal rhoRef = RhoRef->dGet();
477 	if (rhoRef < 0.) {
478 		silent_cerr("StdAirProperties: illegal reference density "
479 			<< rhoRef << std::endl);
480 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
481 	}
482 
483 	if (z < z1) {
484 		/* regione del gradiente (troposfera) */
485 		T = TRef + a * z;
486 		p = PRef * pow(T / TRef, -g0 / (a * R));
487 		rho = rhoRef * pow(T / TRef, -(g0 / (a * R) + 1));
488 
489 	} else {
490 		/* regione a T = const. (stratosfera) */
491 		T = TRef + a * z1;
492 		doublereal p1 = PRef * pow(T / TRef, -g0 / (a * R));
493 		doublereal e = exp(-g0 / (R * T) * (z - z1));
494 		p = p1*e;
495 		doublereal rho1 = rhoRef * pow(T / TRef, -(g0 / (a * R) + 1));
496 		rho = rho1*e;
497 	}
498 
499 	c = sqrt(1.4 * R * T);
500 
501 	return true;
502 }
503 
504 /* StdAirProperties - end */
505 
506 static void
ReadAirstreamData(DataManager * pDM,MBDynParser & HP,TplDriveCaller<Vec3> * & pDC,std::vector<const Gust * > & gust)507 ReadAirstreamData(DataManager *pDM, MBDynParser& HP,
508 		TplDriveCaller<Vec3>*& pDC, std::vector<const Gust*>& gust)
509 {
510 	ASSERT(pDC == 0);
511 
512 	/* Driver multiplo */
513      	pDC = ReadDC3D(pDM, HP);
514 	while (HP.IsKeyWord("gust")) {
515 		/* gust */
516 		gust.insert(gust.end(), ReadGustData(pDM, HP));
517 	}
518 }
519 
520 Elem *
ReadAirProperties(DataManager * pDM,MBDynParser & HP)521 ReadAirProperties(DataManager* pDM, MBDynParser& HP)
522 {
523 	Elem *pEl = NULL;
524 
525 	if (HP.IsKeyWord("std")) {
526 		doublereal PRef(0.);
527 		doublereal rhoRef(0.);
528 		DriveCaller *RhoRef(NULL);
529 		doublereal TRef(0.);
530 		doublereal a(0.);
531 		doublereal R(0.);
532 		doublereal g0(0.);
533 		doublereal z0(0.);
534 		doublereal z1(0.);
535 		doublereal z2(0.);
536 
537 		bool Std = false;
538 
539 		if (HP.IsKeyWord("SI")) {
540 			Std = true;
541 
542 			PRef = 101325.;		/* Pa */
543 			rhoRef = 1.2250;	/* kg/m^3 */
544 			TRef = 288.16;		/* K */
545 			a = -6.5e-3;		/* K/m */
546 			R = 287.;		/* J/kgK */
547 			g0 = 9.81; 		/* m/s^2 */
548 			z1 = 11000.; 		/* m */
549 			z2 = 25000.; 		/* m */
550 
551 		} else if (HP.IsKeyWord("british")) {
552 			Std = true;
553 			PRef = 2116.2; 		/* lb/ft^2 */
554 			rhoRef = 0.002377;	/* slug/ft3 */
555 			TRef = 518.69;		/* R */
556 			a = -3.566e-3;		/* R/ft */
557 			R = 1716;		/* ft lb/slug R */
558 			g0 = 32.17;		/* ft/s^2 */
559 			z1 = 36089;		/* ft */
560 			z2 = 82021;		/* ft */
561 
562 		} else {
563 			PRef = HP.GetReal();
564 			if (PRef <= 0.) {
565 				silent_cerr("illegal reference "
566 					"pressure" << PRef
567 					<< " at line " << HP.GetLineData()
568 					<< std::endl);
569 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
570 			}
571 
572 			RhoRef = HP.GetDriveCaller();
573 			/* FIXME: we need to do runtime checks ... */
574 
575 			TRef = HP.GetReal();
576 			if (TRef <= 0.) {
577 				silent_cerr("illegal reference "
578 					"temperature " << TRef
579 					<< " at line " << HP.GetLineData()
580 					<< std::endl);
581 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
582 			}
583 
584 			a = HP.GetReal();
585 			if (a >= 0.) {
586 				/* FIXME: should we leave this free? */
587 				silent_cerr("illegal temperature gradient "
588 					<< a << " at line " << HP.GetLineData()
589 					<< std::endl);
590 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
591 			}
592 
593 			R = HP.GetReal();
594 			if (R <= 0.) {
595 				silent_cerr("illegal gas constant " << R
596 					<< " at line " << HP.GetLineData()
597 					<< std::endl);
598 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
599 			}
600 
601 			g0 = HP.GetReal();
602 			if (g0 <= 0.) {
603 				silent_cerr("illegal reference "
604 					"gravity acceleration " << g0
605 					<< " at line " << HP.GetLineData()
606 					<< std::endl);
607 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
608 			}
609 
610 			z1 = HP.GetReal();
611 			if (z1 <= 0.) {
612 				silent_cerr("illegal troposphere altitude "
613 					<< z1
614 					<< " at line " << HP.GetLineData()
615 					<< std::endl);
616 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
617 			}
618 
619 			z2 = HP.GetReal();
620 			if (z2 <= z1) {
621 				silent_cerr("illegal stratosphere altitude "
622 					<< z2
623 					<< " at line " << HP.GetLineData()
624 					<< std::endl);
625 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
626 			}
627 		}
628 
629 		if (Std) {
630 			if (HP.IsKeyWord("temperature" "deviation")) {
631 				doublereal T = HP.GetReal();
632 
633 				if (TRef + T <= 0.) {
634 					silent_cerr("illegal "
635 						"temperature deviation " << T
636 						<< " at line "
637 						<< HP.GetLineData()
638 						<< std::endl);
639 					throw ErrGeneric(MBDYN_EXCEPT_ARGS);
640 				}
641 
642 				/*
643 				 * trasformazione isobara applicata
644 				 * all'equazione
645 				 * di stato: rho * R * T = cost
646 				 *
647 				 * rho = rho_0 T_0 / T
648 				 */
649 				rhoRef *= TRef / (TRef + T);
650 				TRef += T;
651 			}
652 
653 			SAFENEWWITHCONSTRUCTOR(RhoRef, ConstDriveCaller,
654 					ConstDriveCaller(rhoRef));
655 		}
656 
657 		if (HP.IsKeyWord("reference" "altitude")) {
658 			z0 = HP.GetReal();
659 		}
660 
661 	     	/* Driver multiplo */
662 	     	TplDriveCaller<Vec3>* pDC = NULL;
663 		std::vector<const Gust *> gust;
664 		ReadAirstreamData(pDM, HP, pDC, gust);
665 	     	flag fOut = pDM->fReadOutput(HP, Elem::AIRPROPERTIES);
666 
667 	     	SAFENEWWITHCONSTRUCTOR(pEl,
668 			StdAirProperties,
669 			StdAirProperties(pDC,
670 				PRef, RhoRef, TRef, a, R, g0, z0, z1, z2,
671 				gust, pDM->pGetRBK(), fOut));
672 	} else {
673 		/* Legacy: density and sound celerity at one altitude;
674 		 * no altitude dependency */
675 
676 		DriveCaller *pRho = HP.GetDriveCaller();
677 
678 	     	doublereal dSS = HP.GetReal();
679 	     	DEBUGLCOUT(MYDEBUG_INPUT, "Sound speed: " << dSS << std::endl);
680 	     	if (dSS <= 0.) {
681 			silent_cerr("illegal null or negative sound speed "
682 				"at line " << HP.GetLineData() << std::endl);
683 
684 			throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
685 	     	}
686 
687 	     	/* Driver multiplo */
688 	     	TplDriveCaller<Vec3>* pDC = NULL;
689 		std::vector<const Gust *> gust;
690 		ReadAirstreamData(pDM, HP, pDC, gust);
691 	     	flag fOut = pDM->fReadOutput(HP, Elem::AIRPROPERTIES);
692 
693 	     	SAFENEWWITHCONSTRUCTOR(pEl,
694 			BasicAirProperties,
695 			BasicAirProperties(pDC, pRho, dSS,
696 				gust, pDM->pGetRBK(), fOut));
697 	}
698 
699 	return pEl;
700 }
701 
702 
703 /* AirPropOwner - begin */
704 
AirPropOwner(void)705 AirPropOwner::AirPropOwner(void)
706 : pAirProperties(NULL)
707 {
708 	NO_OP;
709 }
710 
~AirPropOwner(void)711 AirPropOwner::~AirPropOwner(void)
712 {
713 	NO_OP;
714 }
715 
716 void
PutAirProperties(const AirProperties * pAP)717 AirPropOwner::PutAirProperties(const AirProperties* pAP)
718 {
719 	ASSERT(pAirProperties == NULL);
720 
721 	pAirProperties = pAP;
722 }
723 
724 flag
fGetAirVelocity(Vec3 & Velocity,const Vec3 & X) const725 AirPropOwner::fGetAirVelocity(Vec3& Velocity, const Vec3& X) const
726 {
727 	if (pAirProperties == NULL) {
728 		return 0;
729 	}
730 
731 	Velocity = pAirProperties->GetVelocity(X);
732 	return 1;
733 }
734 
735 doublereal
dGetAirDensity(const Vec3 & X) const736 AirPropOwner::dGetAirDensity(const Vec3& X) const
737 {
738 	return pAirProperties->dGetAirDensity(X);
739 }
740 
741 doublereal
dGetAirPressure(const Vec3 & X) const742 AirPropOwner::dGetAirPressure(const Vec3& X) const
743 {
744 	return pAirProperties->dGetAirPressure(X);
745 }
746 
747 doublereal
dGetAirTemperature(const Vec3 & X) const748 AirPropOwner::dGetAirTemperature(const Vec3& X) const
749 {
750 	return pAirProperties->dGetAirTemperature(X);
751 }
752 
753 doublereal
dGetSoundSpeed(const Vec3 & X) const754 AirPropOwner::dGetSoundSpeed(const Vec3& X) const
755 {
756 	return pAirProperties->dGetSoundSpeed(X);
757 }
758 
759 bool
GetAirProps(const Vec3 & X,doublereal & rho,doublereal & c,doublereal & p,doublereal & T) const760 AirPropOwner::GetAirProps(const Vec3& X, doublereal& rho,
761 		doublereal& c, doublereal& p, doublereal& T) const
762 {
763 	return pAirProperties->GetAirProps(X, rho, c, p, T);
764 }
765 
766 /* AirPropOwner - end */
767 
768 
769 /* AerodynamicElem - begin */
770 
AerodynamicElem(unsigned int uL,const DofOwner * pDO,flag fOut)771 AerodynamicElem::AerodynamicElem(unsigned int uL,
772 	const DofOwner *pDO, flag fOut)
773 : Elem(uL, fOut),
774 ElemWithDofs(uLabel, pDO, fOut)
775 {
776 	NO_OP;
777 }
778 
~AerodynamicElem(void)779 AerodynamicElem::~AerodynamicElem(void)
780 {
781 	NO_OP;
782 }
783 
784 bool
NeedsAirProperties(void) const785 AerodynamicElem::NeedsAirProperties(void) const
786 {
787 	return true;
788 }
789 
790 const InducedVelocity *
pGetInducedVelocity(void) const791 AerodynamicElem::pGetInducedVelocity(void) const
792 {
793 	return 0;
794 }
795 
796 /* AerodynamicElem - end */
797 
798