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