MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
genfm.cc
Go to the documentation of this file.
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
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 
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 
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 
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 
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 
218  = pData->Data[iBeta][iAlpha]*d1Alpha + pData->Data[iBeta][iAlpha + 1]*d2Alpha;
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 
237  const DofOwner *pDO,
238  const StructNode* pN,
239  const Vec3& fTmp, const Mat3x3& RaTmp,
240  doublereal dS, doublereal dL, bool bAlphaFirst,
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 
264 {
265  if (pData != 0) {
266  delete pData;
267  }
268 }
269 
270 /* Tipo dell'elemento (usato per debug ecc.) */
273 {
274  return Elem::AERODYNAMIC;
275 }
276 
277 /* Dimensioni del workspace */
278 void
280 {
281  *piNumRows = 6;
282  *piNumCols = 1;
283 }
284 
285 /* assemblaggio jacobiano */
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
300 {
301  return 0;
302 }
303 
304 /* Dimensioni del workspace */
305 void
307 {
308  *piNumRows = 6;
309  *piNumCols = 1;
310 }
311 
312 /* assemblaggio jacobiano */
315  const VectorHandler& /* XCurr */)
316 {
317  DEBUGCOUTFNAME("GenericAerodynamicForce::InitialAssJac");
318  WorkMat.SetNullMatrix();
319  return WorkMat;
320 }
321 
322 /* Tipo di elemento aerodinamico */
325 {
327 }
328 
329 void
330 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&
338 GenericAerodynamicForce::Restart(std::ostream& out) const
339 {
340  return out;
341 }
342 
343 /* assemblaggio residuo */
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
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
382 {
383  return 2 + 6 + 6;
384 }
385 
386 unsigned int
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 
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:
478  }
479 }
480 
481 /* assemblaggio residuo */
484  const VectorHandler& XCurr)
485 {
486  return WorkVec;
487 }
488 
489 extern Elem *
491  MBDynParser& HP,
492  unsigned int uLabel);
493 
494 /* GenericAerodynamicForce - end */
495 
496 /* GenericAerodynamicData - begin */
497 
498 /* GenericAerodynamicData::GenericAerodynamicCoef - begin */
499 
501 {
502  memset(&dCoef[0], 0, sizeof(dCoef));
503 }
504 
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 
519 {
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 
535 {
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 
550 {
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 
565 {
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 *
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);
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);
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);
622  }
623 
624  if (nBeta <= 0) {
625  silent_cerr("ReadGenericAerodynamicData(" << fname << "): "
626  "invalid number of sideslip angles " << nBeta << " "
627  "from file \"" << fname << "\"" << std::endl);
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);
643  }
644 
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;
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;
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;
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;
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;
795  }
796  }
797 
798  return pData;
799 }
800 
801 Elem *
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);
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);
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);
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);
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);
892 
893  } else {
894  silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
895  "keyword \"file\" or \"reference\" expected "
896  "at line " << HP.GetLineData() << std::endl);
898  }
899 
900  flag fOut = pDM->fReadOutput(HP, Elem::AERODYNAMIC);
901 
902  Elem *pEl = 0;
904  GenericAerodynamicForce(uLabel, pDO,
905  pNode, f, Ra,
906  dRefSurface, dRefLength, bAlphaFirst,
907  pData, fOut));
908 
909  return pEl;
910 }
911 
912 /* GenericAerodynamicData - end */
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
GenericAerodynamicCoef operator-(const GenericAerodynamicCoef &c) const
Definition: genfm.cc:533
virtual VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
Definition: genfm.cc:314
static const doublereal dAlphaMax[]
Definition: genfm.cc:44
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
virtual Elem::Type GetElemType(void) const
Definition: genfm.cc:272
const Vec3 tilde_f
Definition: genfm.h:83
#define M_PI
Definition: gradienttest.cc:67
const Vec3 Zero3(0., 0., 0.)
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
virtual bool GetAirProps(const Vec3 &X, doublereal &rho, doublereal &c, doublereal &p, doublereal &T) const
Definition: aerodyn.cc:760
long int flag
Definition: mbdyn.h:43
void AssVec(SubVectorHandler &WorkVec)
Definition: genfm.cc:49
virtual bool bToBeOutput(void) const
Definition: output.cc:890
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
GenericAerodynamicData * pData
Definition: genfm.h:97
std::vector< doublereal > Alpha
Definition: genfm.h:52
#define DEBUGCOUTFNAME(fname)
Definition: myassert.h:256
virtual void ResizeReset(integer)
Definition: vh.cc:55
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: genfm.cc:345
GenericAerodynamicCoef operator*(const doublereal &d) const
Definition: genfm.cc:549
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
const bool bAlphaFirst
Definition: genfm.h:81
virtual SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: genfm.cc:483
virtual AerodynamicElem::Type GetAerodynamicElemType(void) const
Definition: genfm.cc:324
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
virtual const char * GetFileName(enum Delims Del=DEFAULTDELIM)
Definition: parsinc.cc:673
std::string name
Definition: genfm.h:44
virtual ~GenericAerodynamicForce(void)
Definition: genfm.cc:263
#define NO_OP
Definition: myassert.h:74
const Mat3x3 tilde_Ra
Definition: genfm.h:85
static GenericAerodynamicData * ReadGenericAerodynamicData(const std::string &fname, doublereal dScaleAngle, doublereal dScaleLength, bool bAlphaFirst)
Definition: genfm.cc:582
virtual void InitialWorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: genfm.cc:306
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
std::vector< std::vector< GenericAerodynamicCoef > > Data
Definition: genfm.h:66
doublereal dBeta
Definition: genfm.h:94
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
virtual void Output(OutputHandler &OH) const
Definition: genfm.cc:367
void SetNullMatrix(void)
Definition: submat.h:1159
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
octave_value_list retval
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
GenericAerodynamicCoef operator/(const doublereal &d) const
Definition: genfm.cc:564
std::vector< doublereal > Beta
Definition: genfm.h:53
const doublereal dRefSurface
Definition: genfm.h:79
virtual integer iGetFirstMomentumIndex(void) const =0
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
#define LINE_MAX
Definition: mbdyn.h:163
const StructNode * pNode
Definition: genfm.h:76
const doublereal dRefLength
Definition: genfm.h:80
#define ASSERT(expression)
Definition: colamd.c:977
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
virtual void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
Definition: genfm.cc:330
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
static const doublereal dBetaMax[]
Definition: genfm.cc:45
GenericAerodynamicForce(unsigned int uLabel, const DofOwner *pDO, const StructNode *pN, const Vec3 &fTmp, const Mat3x3 &RaTmp, doublereal dS, doublereal dL, bool bAlphaFirst, GenericAerodynamicData *pData, flag fOut)
Definition: genfm.cc:236
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: genfm.cc:387
static std::stack< cleanup * > c
Definition: cleanup.cc:59
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: genfm.cc:279
const doublereal dS
Definition: beamslider.cc:71
virtual bool GetYesNo(bool &bRet)
Definition: parser.cc:1022
Definition: elem.h:75
Type
Definition: elem.h:91
doublereal dAlpha
Definition: genfm.h:94
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: genfm.cc:287
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
virtual flag fGetAirVelocity(Vec3 &Velocity, const Vec3 &X) const
Definition: aerodyn.cc:725
GradientExpression< UnaryExpr< FuncCos, Expr > > cos(const GradientExpression< Expr > &u)
Definition: gradient.h:2978
Elem * ReadGenericAerodynamicForce(DataManager *pDM, MBDynParser &HP, unsigned int uLabel)
virtual unsigned int iGetInitialNumDof(void) const
Definition: genfm.cc:299
virtual doublereal dGetPrivData(unsigned int i) const
Definition: genfm.cc:442
static doublereal buf[BUFSIZE]
Definition: discctrl.cc:333
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
double doublereal
Definition: colamd.c:52
GenericAerodynamicCoef operator+(const GenericAerodynamicCoef &c) const
Definition: genfm.cc:517
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
virtual std::ostream & Restart(std::ostream &out) const
Definition: genfm.cc:338
virtual unsigned int iGetNumPrivData(void) const
Definition: genfm.cc:381
Mat3x3 R
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056
std::ostream & Aerodynamic(void) const
Definition: output.h:485