MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
rotor.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/rotor.cc,v 1.144 2017/09/26 22:58:54 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 /* Rotore */
33 
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35 
36 #include <limits>
37 #include <cmath>
38 
39 #ifdef USE_MPI
40 #include "mysleep.h"
41 const int mysleeptime = 300;
42 
43 #ifdef MPI_PROFILING
44 extern "C" {
45 #include <mpe.h>
46 #include <stdio.h>
47 }
48 #endif /* MPI_PROFILING */
49 #include "mbcomm.h"
50 #endif /* USE_MPI */
51 
52 #include "rotor.h"
53 #include "dataman.h"
54 
55 static const doublereal dVTipTreshold = 1e-6;
56 
57 /* Rotor - begin */
58 
59 Rotor::Rotor(unsigned int uL, const DofOwner* pDO)
60 : Elem(uL, flag(0)),
61 InducedVelocityElem(uL, pDO, 0, 0, flag(0)),
62 pRotor(0), pGround(0),
63 dOmegaRef(0.), dRadius(0), dVTipRef(0.), dArea(0.),
64 dUMean(0.), dUMeanRef(0.), dUMeanPrev(0.),
65 iMaxIter(0),
66 iCurrIter(0),
67 dTolerance(0),
68 dEta(0),
69 bUMeanRefConverged(false),
70 Weight(), dWeight(0.),
71 dHoverCorrection(1.), dForwardFlightCorrection(1.),
72 RRotTranspose(::Zero3x3), RRot(::Eye3), RRot3(::Zero3),
73 VCraft(::Zero3),
74 dPsi0(0.), dSinAlphad(1.), dCosAlphad(0.),
75 dMu(0.), dLambda(1.), dChi(0.),
76 dVelocity(0.), dOmega(0.),
77 iNumSteps(0)
78 {
79  NO_OP;
80 }
81 
82 Rotor::Rotor(unsigned int uL, const DofOwner* pDO,
83  const StructNode* pC, const Mat3x3& rrot,
84  const StructNode* pR, const StructNode* pG,
85  ResForceSet **ppres,
86  const doublereal& dR,
87  unsigned int iMaxIt, const doublereal& dTol, const doublereal& dE,
88  flag fOut)
89 : Elem(uL, fOut),
90 InducedVelocityElem(uL, pDO, 0, 0, flag(0)),
91 pRotor(0), pGround(0),
92 dOmegaRef(0.), dRadius(0), dVTipRef(0.), dArea(0.),
93 dUMean(0.), dUMeanRef(0.), dUMeanPrev(0.),
94 iMaxIter(0),
95 iCurrIter(0),
96 dTolerance(0),
97 dEta(0),
98 bUMeanRefConverged(false),
99 Weight(), dWeight(0.),
100 dHoverCorrection(1.), dForwardFlightCorrection(1.),
101 RRotTranspose(::Zero3x3), RRot(::Eye3), RRot3(::Zero3),
102 VCraft(::Zero3),
103 dPsi0(0.), dSinAlphad(1.), dCosAlphad(0.),
104 dMu(0.), dLambda(1.), dChi(0.),
105 dVelocity(0.), dOmega(0.),
106 iNumSteps(0)
107 {
108  Init(pC, rrot, pR, pG, ppres, dR, iMaxIt, dTol, dE, fOut);
109 }
110 
111 void
112 Rotor::Init(const StructNode* pC, const Mat3x3& rrot,
113  const StructNode* pR, const StructNode* pG,
114  ResForceSet **ppres,
115  const doublereal& dR,
116  unsigned int iMaxIt, const doublereal& dTol, const doublereal& dE,
117  flag fOut)
118 {
119  InducedVelocityElem::pCraft = pC;
120  InducedVelocityElem::ppRes = ppres;
121 
122  pRotor = pR;
123  pGround = pG;
124  dRadius = dR;
125  iMaxIter = iMaxIt;
126  dTolerance = dTol;
127  dEta = dE;
128  RRot = rrot;
129 
130  SetOutputFlag(fOut);
131 
132  Vec3 R3C(pCraft->GetRCurr()*RRot.GetVec(3));
133  Vec3 R3R(pRotor->GetRCurr().GetVec(3));
134  if (R3C.Dot(R3R) < 1. - std::numeric_limits<doublereal>::epsilon()) {
135  silent_cerr("Rotor(" << GetLabel() << "): "
136  "warning, possible misalignment "
137  "of rotor node StructNode(" << pRotor->GetLabel() << ") "
138  "axis Rr(3) = {" << R3R << "} "
139  "and craft node StructNode(" << pCraft->GetLabel() << ") "
140  "axis Rc*Rh(3) = {" << R3C << "} "
141  << "for Rotor(" << GetLabel() << ")"
142  << std::endl);
143  }
144 }
145 
146 Rotor::~Rotor(void)
147 {
148  NO_OP;
149 }
150 
151 void
152 Rotor::AfterConvergence(const VectorHandler& X,
153  const VectorHandler& XP)
154 {
155  /* non mi ricordo a cosa serve! */
156  iNumSteps++;
157 
158  /* updates the umean at the previous step */
159  dUMeanPrev = dUMean;
160 
161  /*
162  * FIXME: should go in AfterPredict ...
163  * this way there's a one step delay in the weight value
164  */
165  if (Weight.pGetDriveCaller() != 0) {
166  dWeight = Weight.dGet();
167  if (dWeight < 0.) {
168  silent_cout("Rotor(" << GetLabel() << "): "
169  "delay < 0.0; using 0.0" << std::endl);
170  dWeight = 0.;
171 
172  } else if (dWeight > 1.) {
173  silent_cout("Rotor(" << GetLabel() << "): "
174  "delay > 1.0; using 1.0" << std::endl);
175  dWeight = 1.;
176  }
177  }
178 
179  InducedVelocity::AfterConvergence(X, XP);
180 }
181 
182 void
183 Rotor::Output(OutputHandler& OH) const
184 {
185  if (bToBeOutput()) {
186 #ifdef USE_MPI
187  if (is_parallel && IndVelComm.Get_size() > 1) {
188  if (IndVelComm.Get_rank() == 0) {
189  Vec3 TmpF(pTmpVecR);
190  Vec3 TmpM(pTmpVecR+3);
191 
192  OH.Rotors()
193  << std::setw(8) << GetLabel() /* 1 */
194  << " " << RRotTranspose*TmpF /* 2-4 */
195  << " " << RRotTranspose*TmpM /* 5-7 */
196  << " " << dUMean /* 8 */
197  << " " << dVelocity /* 9 */
198  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
199  << " " << dMu /* 11 */
200  << " " << dLambda /* 12 */
201  << " " << dChi /* 13 */
202  << " " << dPsi0 /* 14 */
203  << " " << bUMeanRefConverged /* 15 */
204  << " " << iCurrIter /* 16 */
205  << std::endl;
206 
207  for (int i = 0; ppRes && ppRes[i]; i++) {
208  Vec3 TmpF(pTmpVecR+6+6*i);
209  Vec3 TmpM(pTmpVecR+9+6*i);
210 
211  OH.Rotors()
212  << std::setw(8) << GetLabel()
213  << ":" << ppRes[i]->GetLabel()
214  << " " << TmpF
215  << " " << TmpM
216  << std::endl;
217  }
218  }
219  } else {
220  OH.Rotors()
221  << std::setw(8) << GetLabel() /* 1 */
222  << " " << RRotTranspose*Res.Force() /* 2-4 */
223  << " " << RRotTranspose*Res.Moment() /* 5-7 */
224  << " " << dUMean /* 8 */
225  << " " << dVelocity /* 9 */
226  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
227  << " " << dMu /* 11 */
228  << " " << dLambda /* 12 */
229  << " " << dChi /* 13 */
230  << " " << dPsi0 /* 14 */
231  << " " << bUMeanRefConverged /* 15 */
232  << " " << iCurrIter /* 16 */
233  << std::endl;
234 
235  for (int i = 0; ppRes && ppRes[i]; i++) {
236  OH.Rotors()
237  << std::setw(8) << GetLabel()
238  << ":" << ppRes[i]->GetLabel()
239  << " " << ppRes[i]->pRes->Force()
240  << " " << ppRes[i]->pRes->Moment()
241  << std::endl;
242  }
243  }
244 
245 #else /* !USE_MPI */
246  OH.Rotors()
247  << std::setw(8) << GetLabel() /* 1 */
248  << " " << RRotTranspose*Res.Force() /* 2-4 */
249  << " " << RRotTranspose*Res.Moment() /* 5-7 */
250  << " " << dUMean /* 8 */
251  << " " << dVelocity /* 9 */
252  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
253  << " " << dMu /* 11 */
254  << " " << dLambda /* 12 */
255  << " " << dChi /* 13 */
256  << " " << dPsi0 /* 14 */
257  << " " << bUMeanRefConverged /* 15 */
258  << " " << iCurrIter /* 16 */
259  << std::endl;
260 
261  /* FIXME: check for parallel stuff ... */
262  for (int i = 0; ppRes && ppRes[i]; i++) {
263  OH.Rotors()
264  << std::setw(8) << GetLabel()
265  << ":" << ppRes[i]->GetLabel()
266  << " " << ppRes[i]->pRes->Force()
267  << " " << ppRes[i]->pRes->Moment()
268  << std::endl;
269  }
270 #endif /* !USE_MPI */
271  }
272 }
273 
274 /* Calcola la posizione azimuthale di un punto generico.
275  * X e' il punto di cui e' chiesta la posizione azimuthale,
276  * nel sistema assoluto.
277  * Gli viene sottratta la posizione del rotore (corpo attorno al quale avviene
278  * la rotazione). Quindi il vettore risultante viene trasformato
279  * nel riferimento del mozzo.
280  * Si trascura l'eventuale componente fuori del piano, dalle altre due si
281  * ricava l'angolo di rotazione relativa, a cui viene sommato l'angolo
282  * che il corpo rotore forma con la direzione del vento relativo dPsi0,
283  * calcolata in precedenza.
284  */
285 doublereal
286 Rotor::dGetPsi(const Vec3& X) const
287 {
288  doublereal dr, dp;
289  GetPos(X, dr, dp);
290  return dp;
291 }
292 
293 /* Calcola la distanza di un punto dall'asse di rotazione in coordinate
294  * adimensionali */
295 doublereal
296 Rotor::dGetPos(const Vec3& X) const
297 {
298  doublereal dr, dp;
299  GetPos(X, dr, dp);
300  return dr;
301 }
302 
303 void
304 Rotor::GetPos(const Vec3& X, doublereal& dr, doublereal& dp) const
305 {
306  Vec3 XRel(RRotTranspose*(X - Res.Pole()));
307 
308  doublereal d1 = XRel(1);
309  doublereal d2 = XRel(2);
310 
311  doublereal d = sqrt(d1*d1 + d2*d2);
312 
313  ASSERT(dRadius > std::numeric_limits<doublereal>::epsilon());
314  dr = d/dRadius;
315 
316  dp = atan2(d2, d1) - dPsi0;
317 }
318 
319 /* Calcola vari parametri geometrici
320  * A partire dai corpi che identificano il velivolo ed il rotore
321  */
322 void
323 Rotor::InitParam(bool bComputeMeanInducedVelocity)
324 {
325  ASSERT(pCraft != 0);
326  ASSERT(pRotor != 0);
327 
328  /* Trasposta della matrice di rotazione del rotore */
329  RRotTranspose = pCraft->GetRCurr()*RRot;
330  RRot3 = RRotTranspose.GetVec(3);
331  RRotTranspose = RRotTranspose.Transpose();
332 
333  /* Posizione del rotore */
334  Res.PutPole(pRotor->GetXCurr());
335 
336  /* Velocita' angolare del rotore */
337  dOmega = (pRotor->GetWCurr() - pCraft->GetWCurr()).Norm();
338 
339  /* Velocita' di traslazione del velivolo */
340  VCraft = -pRotor->GetVCurr();
341  Vec3 VTmp(Zero3);
342  if (fGetAirVelocity(VTmp, pRotor->GetXCurr())) {
343  VCraft += VTmp;
344  }
345 
346  /* Velocita' nel sistema del velivolo (del disco?) decomposta */
347  VTmp = RRotTranspose*VCraft;
348  doublereal dV1 = VTmp(1);
349  doublereal dV2 = VTmp(2);
350  doublereal dV3 = VTmp(3);
351  doublereal dVV = dV1*dV1 + dV2*dV2;
352  doublereal dV = sqrt(dVV);
353 
354  /* Angolo di azimuth 0 del rotore */
355  dPsi0 = atan2(dV2, dV1);
356 
357  /* Angolo di influsso */
358  dVelocity = sqrt(dV3*dV3 + dVV);
359  if (dVelocity > dVTipTreshold*dVTipRef) {
360  dSinAlphad = -dV3/dVelocity;
361  dCosAlphad = dV/dVelocity;
362 
363  } else {
364  dSinAlphad = 1.;
365  dCosAlphad = 0.;
366  }
367 
368  if (!bComputeMeanInducedVelocity) {
369  return;
370  }
371 
372  // bail out when density is negligible
373  doublereal dRho = dGetAirDensity(GetXCurr());
374  if (dRho <= std::numeric_limits<doublereal>::epsilon()) {
375  return;
376  }
377 
378  // Thrust in rotor reference frame
379  doublereal dT = RRot3*Res.Force();
380 
381  /* Parametri di influsso (usano il valore di dUMean al passo precedente) */
382  doublereal dVTip = 0.;
383  dMu = 0.;
384  dLambda = 0.;
385  dVTip = dOmega*dRadius;
386 
387  if (dVTip > dVTipTreshold*dVTipRef) {
388  dMu = (dVelocity*dCosAlphad)/dVTip;
389  }
390 
391  /* NOTE: dUMeanRef starts at the value it had previously */
392  if (std::abs(dUMeanRef) < std::numeric_limits<doublereal>::epsilon()
393  && std::abs(dT) > std::numeric_limits<doublereal>::epsilon())
394  {
395  // first guess: start with the value in hover
396  dUMeanRef = copysign(std::sqrt(std::abs(dT)/(2*dRho*dArea)), dT);
397  }
398 
399  /* Ground effect */
400  doublereal dGE = 1.;
401  if (pGround) {
402  doublereal p = pGround->GetRCurr().GetVec(3)*(pRotor->GetXCurr() - pGround->GetXCurr());
403  doublereal z = p/(dRadius/4.);
404 
405  if (z < .25) {
406  if (z <= 0.) {
407  silent_cerr("warning, illegal negative "
408  "normalized altitude "
409  "z=" << z << std::endl);
410  }
411 
412  z = .25;
413  }
414 
415  /*
416  * According to I.C. Cheeseman & N.E. Bennett,
417  * "The Effect of Ground on a Helicopter Rotor
418  * in Forward Flight", NASA TR-3021, 1955:
419  *
420  * U = Uref * ( 1 - ( R / ( 4 * z ) )^2 )
421  *
422  * We need to make R / ( 4 * z ) <= 1, so
423  * we must enforce z >= R / 4.
424  */
425  dGE -= 1./(z*z);
426  }
427 
428  bUMeanRefConverged = false;
429  if (dVTip > dVTipTreshold*dVTipRef) {
430  doublereal dCt = dT/(dRho*dArea*dVTip*dVTip);
431  doublereal dLambda0 = dVelocity*dSinAlphad/dVTip;
432  doublereal dLambdaInd = dUMeanRef/dVTip;
433 
434  bool bRetrying = false;
435 retry:;
436 
437  for (iCurrIter = 0; iCurrIter < iMaxIter; iCurrIter++) {
438  dLambda = dLambda0 + dLambdaInd;
439 
440  doublereal dRef0 = dMu*dMu + dLambda*dLambda;
441  doublereal dRef1 = 2.*sqrt(dRef0);
442  doublereal dRef2 = dRef1*dRef0;
443  doublereal dF = 0.;
444  if (dRef1 > std::numeric_limits<doublereal>::epsilon()) {
445  dF = dLambdaInd - dCt/dRef1;
446  doublereal dFPrime = 1. + (dCt/dRef2)*dLambda;
447  doublereal dDelta = dF/dFPrime;
448  dLambdaInd -= dEta*dDelta;
449  }
450 
451 #if 0
452  std::cerr
453  << " iter=" << iCurrIter
454  << " lambda=" << dLambda
455  << " dRef1=" << dRef1
456  << " dF=" << dF
457  << " Uref=" << dUMeanRef
458  << std::endl;
459 #endif
460 
461  if (std::abs(dF) <= dTolerance) {
462  bUMeanRefConverged = true;
463  break;
464  }
465  }
466 
467  // NOTE: the induced velocity must have the same sign
468  // of the thrust coefficient
469  if (dLambdaInd*dCt < 0.) {
470  if (!bRetrying) {
471  bRetrying = true;
472 
473  // first guess: revert the sign of induced velocity
474  dLambdaInd = -dLambdaInd;
475  goto retry;
476  }
477 
478  silent_cerr("Rotor(" << GetLabel() << "): "
479  "induced velocity and thrust signs differ "
480  "(lambda_u=" << dLambdaInd << ", Ct=" << dCt << ")"
481  << std::endl);
482  }
483 
484  dLambda = dLambda0 + dLambdaInd;
485 
486  // this is updated to serve as starting point for next iteration
487  dUMeanRef = dLambdaInd*dVTip;
488 
489  /* if no convergence, simply accept the current value
490  * very forgiving choice, though */
491 #if 0
492  if (iCurrIter == iMaxIter) {
493  silent_cerr("unable to compute mean induced velocity "
494  "for Rotor(" << GetLabel() << ")" << std::endl);
495  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
496  }
497 #endif
498  }
499 
500  //if (dMu == 0. && dLambda == 0.) {
501  if (std::abs(dMu) < std::numeric_limits<doublereal>::epsilon() ) {
502  dChi = 0.;
503 
504  } else {
505  dChi = atan2(dMu, dLambda);
506  }
507 
508  /*
509  * From Claudio Monteggia:
510  *
511  * Ct
512  * Um = -------------------------------------
513  * sqrt( lambda^2 / KH^4 + mu^2 / KF^2)
514  */
515  if (dVTip > dVTipTreshold*dVTipRef) {
516  doublereal dMuTmp = dMu/dForwardFlightCorrection;
517  doublereal dLambdaTmp = dLambda/(dHoverCorrection*dHoverCorrection);
518  doublereal dRef = 2*sqrt(dMuTmp*dMuTmp + dLambdaTmp*dLambdaTmp);
519  doublereal dCt = dT/(dRho*dArea*dVTip*dVTip);
520  if (dRef > std::abs(dCt)) {
521  // NOTE: an inflow velocity larger than VTip
522  // makes little sense
523  dUMean = (1. - dWeight)*dGE*dVTip*dCt/dRef + dWeight*dUMeanPrev;
524 
525  } else {
526  dUMean = dUMeanPrev;
527  }
528 
529  } else {
530  dUMean = dUMeanPrev;
531  }
532 }
533 
534 std::ostream&
535 Rotor::Restart(std::ostream& out) const
536 {
537  return out << " rotor: " << GetLabel() << ", "
538  << pCraft->GetLabel() << ", " << pRotor->GetLabel()
539  << ", induced velocity: ";
540 }
541 
542 /* Rotor - end */
543 
544 
545 /* NoRotor - begin */
546 
547 NoRotor::NoRotor(unsigned int uLabel,
548  const DofOwner* pDO)
549 : Elem(uLabel, flag(0)),
550 Rotor(uLabel, pDO)
551 {
552  NO_OP;
553 }
554 
555 NoRotor::NoRotor(unsigned int uLabel,
556  const DofOwner* pDO,
557  const StructNode* pCraft,
558  const Mat3x3& rrot,
559  const StructNode* pRotor,
560  ResForceSet **ppres,
561  const doublereal& dR,
562  flag fOut)
563 : Elem(uLabel, flag(0)),
564 Rotor(uLabel, pDO)
565 {
566  Init(pCraft, rrot, pRotor, ppres, dR, fOut);
567 }
568 
569 void
570 NoRotor::Init(const StructNode* pCraft,
571  const Mat3x3& rrot,
572  const StructNode* pRotor,
573  ResForceSet **ppres,
574  const doublereal& dR,
575  flag fOut)
576 {
577  Rotor::Init(pCraft, rrot, pRotor, 0, ppres, dR, 0, 0., 0., fOut);
578 
579 #ifdef USE_MPI
580  if (is_parallel && bToBeOutput()) {
581  SAFENEWARR(pBlockLenght, int, 3);
582  SAFENEWARR(pDispl, MPI::Aint, 3);
583  for (int i = 0; i < 3; i++) {
584  pBlockLenght[i] = 1;
585  }
586  for (int i = 0; i < 3; i++) {
587  pDispl[i] = MPI::Get_address(&(Res.Pole().pGetVec()[i]));
588  }
589  SAFENEWWITHCONSTRUCTOR(pIndVelDataType, MPI::Datatype,
590  MPI::Datatype(MPI::DOUBLE.Create_hindexed(3, pBlockLenght, pDispl)));
591  pIndVelDataType->Commit();
592  }
593 #endif /* USE_MPI */
594 }
595 
596 NoRotor::~NoRotor(void)
597 {
598 #ifdef USE_MPI
599  SAFEDELETEARR(pBlockLenght);
600  SAFEDELETEARR(pDispl);
601  SAFEDELETE(pIndVelDataType);
602 #endif /* USE_MPI */
603 }
604 
605 /* assemblaggio residuo */
606 SubVectorHandler&
607 NoRotor::AssRes(SubVectorHandler& WorkVec,
608  doublereal /* dCoef */ ,
609  const VectorHandler& /* XCurr */ ,
610  const VectorHandler& /* XPrimeCurr */ )
611 {
612  DEBUGCOUT("Entering NoRotor::AssRes()" << std::endl);
613 
614  flag out = fToBeOutput();
615 
616 #ifdef USE_MPI
617  if (out) {
618  ExchangeLoads(out);
619  }
620 
621  if (!is_parallel || IndVelComm.Get_rank() == 0)
622 #endif /* USE_MPI */
623  {
624  if (out) {
625  /* Calcola parametri vari */
626  Rotor::InitParam(false);
627 
628 #ifdef USE_MPI
629  ExchangeVelocity();
630 #endif /* USE_MPI */
631  }
632  }
633 
634  ResetForce();
635  WorkVec.Resize(0);
636 
637 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
638  Done();
639 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
640 
641  return WorkVec;
642 }
643 
644 std::ostream&
645 NoRotor::Restart(std::ostream& out) const
646 {
647  return Rotor::Restart(out) << "no;" << std::endl;
648 }
649 
650 // Somma alla trazione il contributo di forza di un elemento generico
651 void
652 NoRotor::AddForce(const Elem *pEl, const StructNode *pNode,
653  const Vec3& F, const Vec3& M, const Vec3& X)
654 {
655  // Non gli serve in quanto non calcola velocita' indotta.
656  // Solo se deve fare l'output lo calcola
657 #ifdef USE_MPI
658  if (ReqV != MPI::REQUEST_NULL) {
659  while (!ReqV.Test()) {
660  MYSLEEP(mysleeptime);
661  }
662  }
663 #endif // USE_MPI
664 
665 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
666  pthread_mutex_lock(&forces_mutex);
667  Wait();
668 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
669 
670  if (bToBeOutput()) {
671  Res.AddForces(F, M, X);
672  InducedVelocity::AddForce(pEl, pNode, F, M, X);
673  }
674 
675 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
676  pthread_mutex_unlock(&forces_mutex);
677 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
678 }
679 
680 /* Restituisce ad un elemento la velocita' indotta in base alla posizione
681  * azimuthale */
682 Vec3
683 NoRotor::GetInducedVelocity(Elem::Type type,
684  unsigned uLabel, unsigned uPnt, const Vec3& X) const
685 {
686  return Zero3;
687 }
688 
689 /* NoRotor - end */
690 
691 
692 /* UniformRotor - begin */
693 
694 UniformRotor::UniformRotor(unsigned int uLabel, const DofOwner* pDO)
695 : Elem(uLabel, flag(0)),
696 Rotor(uLabel, pDO)
697 {
698  NO_OP;
699 }
700 
701 UniformRotor::UniformRotor(unsigned int uLabel,
702  const DofOwner* pDO,
703  const StructNode* pCraft,
704  const Mat3x3& rrot,
705  const StructNode* pRotor,
706  const StructNode* pGround,
707  ResForceSet **ppres,
708  const doublereal& dOR,
709  const doublereal& dR,
710  DriveCaller *pdW,
711  unsigned int iMaxIt,
712  const doublereal& dTol,
713  const doublereal& dE,
714  const doublereal& dCH,
715  const doublereal& dCFF,
716  flag fOut)
717 : Elem(uLabel, flag(0)),
718 Rotor(uLabel, pDO)
719 {
720  Init(pCraft, rrot, pRotor, pGround, ppres, dOR, dR,
721  pdW, iMaxIt, dTol, dE, dCH, dCFF, fOut);
722 }
723 
724 void
725 UniformRotor::Init(const StructNode* pCraft,
726  const Mat3x3& rrot,
727  const StructNode* pRotor,
728  const StructNode* pGround,
729  ResForceSet **ppres,
730  const doublereal& dOR,
731  const doublereal& dR,
732  DriveCaller *pdW,
733  unsigned int iMaxIt,
734  const doublereal& dTol,
735  const doublereal& dE,
736  const doublereal& dCH,
737  const doublereal& dCFF,
738  flag fOut)
739 {
740  ASSERT(dOR > 0.);
741  ASSERT(dR > 0.);
742  ASSERT(pdW != 0);
743 
744  Rotor::Init(pCraft, rrot, pRotor, pGround, ppres, dR, iMaxIt, dTol, dE, fOut);
745 
746  dOmegaRef = dOR;
747  dVTipRef = dOmegaRef*dRadius;
748  dArea = M_PI*dRadius*dRadius;
749  Weight.Set(pdW);
750  dHoverCorrection = dCH;
751  dForwardFlightCorrection = dCFF;
752 
753 #ifdef USE_MPI
754  if (is_parallel) {
755  SAFENEWARR(pBlockLenght, int, 7);
756  SAFENEWARR(pDispl, MPI::Aint, 7);
757  for (int i = 0; i < 7; i++) {
758  pBlockLenght[i] = 1;
759  }
760  for (int i = 0; i < 3; i++) {
761  pDispl[i] = MPI::Get_address(&(RRot3.pGetVec()[i]));
762  }
763  pDispl[3] = MPI::Get_address(&dUMeanPrev);
764  for (int i = 4; i <= 6; i++) {
765  pDispl[i] = MPI::Get_address(&(Res.Pole().pGetVec()[i-4]));
766  }
767  SAFENEWWITHCONSTRUCTOR(pIndVelDataType, MPI::Datatype,
768  MPI::Datatype(MPI::DOUBLE.Create_hindexed(7, pBlockLenght, pDispl)));
769  pIndVelDataType->Commit();
770  }
771 #endif /* USE_MPI */
772 }
773 
774 UniformRotor::~UniformRotor(void)
775 {
776 #ifdef USE_MPI
777  SAFEDELETEARR(pBlockLenght);
778  SAFEDELETEARR(pDispl);
779  SAFEDELETE(pIndVelDataType);
780 #endif /* USE_MPI */
781 }
782 
783 /* assemblaggio residuo */
784 SubVectorHandler&
785 UniformRotor::AssRes(SubVectorHandler& WorkVec,
786  doublereal /* dCoef */ ,
787  const VectorHandler& /* XCurr */ ,
788  const VectorHandler& /* XPrimeCurr */ )
789 {
790  DEBUGCOUT("Entering UniformRotor::AssRes()" << std::endl);
791 
792 #ifdef USE_MPI
793  ExchangeLoads(fToBeOutput());
794  if (!is_parallel || IndVelComm.Get_rank() == 0)
795 #endif /* USE_MPI */
796  {
797  /* Calcola parametri vari */
798  Rotor::InitParam();
799 
800 #ifdef DEBUG
801  // Prova:
802 #if 0
803  Vec3 XTmp(2.,2.,0.);
804  doublereal dPsiTmp = dGetPsi(XTmp);
805  doublereal dXTmp = dGetPos(XTmp);
806  std::cout
807  << "X rotore: " << pRotor->GetXCurr() << std::endl
808  << "V rotore: " << VCraft << std::endl
809  << "X punto: " << XTmp << std::endl
810  << "Omega: " << dOmega << std::endl
811  << "Velocita': " << dVelocity << std::endl
812  << "Psi0: " << dPsi0 << std::endl
813  << "Psi punto: " << dPsiTmp << std::endl
814  << "Raggio: " << dRadius << std::endl
815  << "r punto: " << dXTmp << std::endl
816  << "mu: " << dMu << std::endl
817  << "lambda: " << dLambda << std::endl
818  << "cos(ad): " << dCosAlphad << std::endl
819  << "sin(ad): " << dSinAlphad << std::endl
820  << "UMean: " << dUMean << std::endl;
821 #endif
822 #endif /* DEBUG */
823  }
824 
825 #ifdef USE_MPI
826  ExchangeVelocity();
827 #endif /* USE_MPI */
828 
829  ResetForce();
830 
831  /* Non tocca il residuo */
832  WorkVec.Resize(0);
833 
834 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
835  Done();
836 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
837 
838  return WorkVec;
839 }
840 
841 std::ostream&
842 UniformRotor::Restart(std::ostream& out) const
843 {
844  return Rotor::Restart(out) << "uniform, " << dRadius << ", ",
845  Weight.pGetDriveCaller()->Restart(out)
846  << ", correction, " << dHoverCorrection
847  << ", " << dForwardFlightCorrection << ';' << std::endl;
848 }
849 
850 /* Somma alla trazione il contributo di forza di un elemento generico */
851 void
852 UniformRotor::AddForce(const Elem *pEl, const StructNode *pNode,
853  const Vec3& F, const Vec3& M, const Vec3& X)
854 {
855 #ifdef USE_MPI
856  if (ReqV != MPI::REQUEST_NULL) {
857  while (!ReqV.Test()) {
858  MYSLEEP(mysleeptime);
859  }
860  }
861 #endif /* USE_MPI */
862 
863 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
864  pthread_mutex_lock(&forces_mutex);
865  Wait();
866 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
867 
868  /* Solo se deve fare l'output calcola anche il momento */
869  if (bToBeOutput()) {
870  Res.AddForces(F, M, X);
871  InducedVelocity::AddForce(pEl, pNode, F, M, X);
872  } else {
873  Res.AddForce(F);
874  }
875 
876 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
877  pthread_mutex_unlock(&forces_mutex);
878 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
879 }
880 
881 /* Restituisce ad un elemento la velocita' indotta in base alla posizione
882  * azimuthale */
883 Vec3
884 UniformRotor::GetInducedVelocity(Elem::Type type,
885  unsigned uLabel, unsigned uPnt, const Vec3& X) const
886 {
887 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
888  Wait();
889 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
890 
891  return RRot3*dUMeanPrev;
892 };
893 
894 UniformRotor2::UniformRotor2(unsigned int uLabel, const DofOwner* pDO)
895 : Elem(uLabel, flag(0)),
896 UniformRotor(uLabel, pDO)
897 {
898  NO_OP;
899 }
900 
901 UniformRotor2::UniformRotor2(unsigned int uLabel,
902  const DofOwner* pDO,
903  const StructNode* pCraft,
904  const Mat3x3& rrot,
905  const StructNode* pRotor,
906  const StructNode* pGround,
907  ResForceSet **ppres,
908  const doublereal& dOR,
909  const doublereal& dR,
910  DriveCaller *pdW,
911  unsigned int iMaxIt,
912  const doublereal& dTol,
913  const doublereal& dE,
914  const doublereal& dCH,
915  const doublereal& dCFF,
916  flag fOut)
917 : Elem(uLabel, fOut),
918 UniformRotor(uLabel, pDO, pCraft, rrot, pRotor, pGround, ppres, dOR, dR, pdW, iMaxIt, dTol, dE, dCH, dCFF, fOut)
919 {
920  NO_OP;
921 }
922 
923 UniformRotor2::~UniformRotor2(void)
924 {
925  NO_OP;
926 }
927 
928 bool
929 UniformRotor2::bSectionalForces(void) const
930 {
931  return true;
932 }
933 
934 /* Somma alla trazione il contributo di forza di un elemento generico */
935 void
936 UniformRotor2::AddSectionalForce(Elem::Type type,
937  const Elem *pEl, unsigned uPnt,
938  const Vec3& F, const Vec3& M, doublereal dW,
939  const Vec3& X, const Mat3x3& R,
940  const Vec3& V, const Vec3& W)
941 {
942 #ifdef USE_MPI
943  if (ReqV != MPI::REQUEST_NULL) {
944  while (!ReqV.Test()) {
945  MYSLEEP(mysleeptime);
946  }
947  }
948 #endif /* USE_MPI */
949 
950 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
951  pthread_mutex_lock(&forces_mutex);
952  Wait();
953 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
954 
955  /* Solo se deve fare l'output calcola anche il momento */
956  if (bToBeOutput()) {
957  Vec3 FTmp(F*dW);
958  Vec3 MTmp(M*dW);
959  Res.AddForces(FTmp, MTmp, X);
960  InducedVelocity::AddForce(pEl, 0, FTmp, MTmp, X);
961 
962  } else {
963  Res.AddForce(F*dW);
964  }
965 
966 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
967  pthread_mutex_unlock(&forces_mutex);
968 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
969 }
970 
971 /* UniformRotor - end */
972 
973 
974 /* GlauertRotor - begin */
975 
976 GlauertRotor::GlauertRotor(unsigned int uLabel, const DofOwner* pDO)
977 : Elem(uLabel, flag(0)),
978 Rotor(uLabel, pDO),
979 gtype(GlauertRotor::UNKNOWN)
980 {
981  NO_OP;
982 }
983 
984 GlauertRotor::GlauertRotor(unsigned int uLabel,
985  const DofOwner* pDO,
986  const StructNode* pCraft,
987  const Mat3x3& rrot,
988  const StructNode* pRotor,
989  const StructNode* pGround,
990  ResForceSet **ppres,
991  const doublereal& dOR,
992  const doublereal& dR,
993  DriveCaller *pdW,
994  unsigned int iMaxIt,
995  const doublereal& dTol,
996  const doublereal& dE,
997  const doublereal& dCH,
998  const doublereal& dCFF,
999  GlauertRotor::Type gtype,
1000  flag fOut)
1001 : Elem(uLabel, flag(0)),
1002 Rotor(uLabel, pDO),
1003 gtype(gtype)
1004 {
1005  Init(pCraft, rrot, pRotor, pGround, ppres, dOR, dR,
1006  pdW, iMaxIt, dTol, dE, dCH, dCFF, fOut);
1007 }
1008 
1009 void
1010 GlauertRotor::Init(const StructNode* pCraft,
1011  const Mat3x3& rrot,
1012  const StructNode* pRotor,
1013  const StructNode* pGround,
1014  ResForceSet **ppres,
1015  const doublereal& dOR,
1016  const doublereal& dR,
1017  DriveCaller *pdW,
1018  unsigned int iMaxIt,
1019  const doublereal& dTol,
1020  const doublereal& dE,
1021  const doublereal& dCH,
1022  const doublereal& dCFF,
1023  flag fOut)
1024 {
1025  ASSERT(dOR > 0.);
1026  ASSERT(dR > 0.);
1027  ASSERT(pdW != 0);
1028 
1029  Rotor::Init(pCraft, rrot, pRotor, pGround, ppres, dR, iMaxIt, dTol, dE, fOut);
1030 
1031  dOmegaRef = dOR;
1032  dVTipRef = dOmegaRef*dRadius;
1033  dArea = M_PI*dRadius*dRadius;
1034  Weight.Set(pdW);
1035  dHoverCorrection = dCH;
1036  dForwardFlightCorrection = dCFF;
1037 
1038 #ifdef USE_MPI
1039  if (is_parallel) {
1040  SAFENEWARR(pBlockLenght, int, 20);
1041  SAFENEWARR(pDispl, MPI::Aint, 20);
1042  for (int i = 0; i < 20; i++) {
1043  pBlockLenght[i] = 1;
1044  }
1045  for (int i = 0; i < 3; i++) {
1046  pDispl[i] = MPI::Get_address(&(RRot3.pGetVec()[i]));
1047  }
1048  pDispl[3] = MPI::Get_address(&dUMeanPrev);
1049  pDispl[4] = MPI::Get_address(&dLambda);
1050  pDispl[5] = MPI::Get_address(&dMu);
1051  pDispl[6] = MPI::Get_address(&dChi);
1052  pDispl[7] = MPI::Get_address(&dPsi0);
1053  for (int i = 8; i <= 10; i++) {
1054  pDispl[i] = MPI::Get_address(&(Res.Pole().pGetVec()[i-8]));
1055  }
1056  for (int i = 11; i < 20; i++) {
1057  pDispl[i] = MPI::Get_address(&(RRotTranspose.pGetMat()[i-11]));
1058  }
1059  SAFENEWWITHCONSTRUCTOR(pIndVelDataType, MPI::Datatype,
1060  MPI::Datatype(MPI::DOUBLE.Create_hindexed(20, pBlockLenght, pDispl)));
1061  pIndVelDataType->Commit();
1062  }
1063 #endif /* USE_MPI */
1064 }
1065 
1066 
1067 GlauertRotor::~GlauertRotor(void)
1068 {
1069 #ifdef USE_MPI
1070  SAFEDELETEARR(pBlockLenght);
1071  SAFEDELETEARR(pDispl);
1072  SAFEDELETE(pIndVelDataType);
1073 #endif /* USE_MPI */
1074 }
1075 
1076 
1077 /* assemblaggio residuo */
1078 SubVectorHandler&
1079 GlauertRotor::AssRes(SubVectorHandler& WorkVec,
1080  doublereal /* dCoef */ ,
1081  const VectorHandler& /* XCurr */ ,
1082  const VectorHandler& /* XPrimeCurr */ )
1083 {
1084  DEBUGCOUT("Entering GlauertRotor::AssRes()" << std::endl);
1085 
1086 #ifdef USE_MPI
1087  ExchangeLoads(fToBeOutput());
1088  if (!is_parallel || IndVelComm.Get_rank() == 0)
1089 #endif /* USE_MPI */
1090  {
1091  /* Calcola parametri vari */
1092  Rotor::InitParam();
1093  }
1094 
1095 #ifdef USE_MPI
1096  ExchangeVelocity();
1097 #endif /* USE_MPI */
1098 
1099  ResetForce();
1100 
1101  /* Non tocca il residuo */
1102  WorkVec.Resize(0);
1103 
1104 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1105  Done();
1106 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
1107 
1108  return WorkVec;
1109 }
1110 
1111 std::ostream&
1112 GlauertRotor::Restart(std::ostream& out) const
1113 {
1114  return Rotor::Restart(out) << "Glauert, " << dRadius << ", ",
1115  Weight.pGetDriveCaller()->Restart(out)
1116  << ", correction, " << dHoverCorrection
1117  << ", " << dForwardFlightCorrection << ';' << std::endl;
1118 }
1119 
1120 
1121 /* Somma alla trazione il contributo di forza di un elemento generico */
1122 void
1123 GlauertRotor::AddForce(const Elem *pEl, const StructNode *pNode,
1124  const Vec3& F, const Vec3& M, const Vec3& X)
1125 {
1126 #ifdef USE_MPI
1127  if (ReqV != MPI::REQUEST_NULL) {
1128  while (!ReqV.Test()) {
1129  MYSLEEP(mysleeptime);
1130  }
1131  }
1132 #endif /* USE_MPI */
1133 
1134 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1135  pthread_mutex_lock(&forces_mutex);
1136  Wait();
1137 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
1138 
1139  /* Solo se deve fare l'output calcola anche il momento */
1140  if (bToBeOutput()) {
1141  Res.AddForces(F, M, X);
1142  InducedVelocity::AddForce(pEl, pNode, F, M, X);
1143  } else {
1144  Res.AddForce(F);
1145  }
1146 
1147 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1148  pthread_mutex_unlock(&forces_mutex);
1149 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
1150 }
1151 
1152 
1153 /*
1154  * Restituisce ad un elemento la velocita' indotta in base alla posizione
1155  * azimuthale
1156  */
1157 Vec3
1158 GlauertRotor::GetInducedVelocity(Elem::Type type,
1159  unsigned uLabel, unsigned uPnt, const Vec3& X) const
1160 {
1161  if (dUMeanPrev == 0.) {
1162  return Zero3;
1163  }
1164 
1165 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1166  Wait();
1167 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
1168 
1169  if (std::abs(dLambda) < 1.e-9) {
1170  return RRot3*dUMeanPrev;
1171  }
1172 
1173  doublereal dr, dp;
1174  GetPos(X, dr, dp);
1175 
1176  // recall that tan(dChi) = dMu/dLambda
1177 
1178  doublereal k1, k2 = 0.;
1179  switch (gtype) {
1180  case GLAUERT:
1181  // NASA-TM-102219: attributed to Drees et al.
1182  k1 = 4./3.*(1. - 1.8*dMu*dMu)*tan(dChi/2.);
1183  break;
1184 
1185 #if 0
1186  case WHEATLEY:
1187  // NASA-TM-102219
1188  k1 = 0.5;
1189  break;
1190 #endif
1191 
1192  case COLEMAN_ET_AL:
1193  // NASA-TM-102219
1194  // Gordon J. Leishman, "Principles of Helicopter Aerodynamics", 2nd Ed., 2006
1195  // Wayne Johnson, "Rotorcraft Aeromechanics", 2013
1196  k1 = tan(dChi/2.);
1197  break;
1198 
1199  case DREES_1:
1200  case DREES_2:
1201  // Gordon J. Leishman, "Principles of Helicopter Aerodynamics", 2nd Ed., 2006
1202  // Wayne Johnson, "Rotorcraft Aeromechanics", 2013
1203  // k1 = 4./3.*(1 - cos(dChi) - 1.8*dMu*dMu)/sin(dChi); // risk of division by zero...
1204  k1 = 4./3.*(tan(dChi/2.) - 1.8*dLambda*dMu/cos(dChi));
1205  k2 = -2.*dMu;
1206  break;
1207 
1208  case PAYNE:
1209  // NASA-TM-102219
1210  // Gordon J. Leishman, "Principles of Helicopter Aerodynamics", 2nd Ed., 2006
1211  // k1 = 4./3.*(dMu/dLambda/(1.2 + dMu/dLambda));
1212  // reduce risk of division by zero
1213  k1 = 4./3.*dMu/(1.2*dLambda + dMu);
1214  break;
1215 
1216  case WHITE_AND_BLAKE:
1217  // NASA-TM-102219
1218  // Gordon J. Leishman, "Principles of Helicopter Aerodynamics", 2nd Ed., 2006
1219  // Wayne Johnson, "Rotorcraft Aeromechanics", 2013
1220  k1 = sqrt(2.)*sin(dChi);
1221  break;
1222 
1223  case PITT_AND_PETERS:
1224  // NASA-TM-102219
1225  // Gordon J. Leishman, "Principles of Helicopter Aerodynamics", 2nd Ed., 2006: "23" instead of "32"...
1226  k1 = 15.*M_PI/32.*tan(dChi/2.);
1227  break;
1228 
1229  case HOWLETT:
1230  // NASA-TM-102219
1231  // Gordon J. Leishman, "Principles of Helicopter Aerodynamics", 2nd Ed., 2006
1232  k1 = pow(sin(dChi), 2);
1233  break;
1234 
1235 #if 0
1236  case DREES_2: {
1237  // Jianhua Zhang, "Active-Passive Hybrid Optimization of Rotor Blades With Trailing Edge Flaps", PhD Thesis, PSU, 2001
1238  // in the end, it is identical to Drees' as of Wayne Johnson and Gordon J. Leishman
1239  // FIXME: what if dMu ~ 0?
1240  doublereal dLdM = dLambda/dMu;
1241  // k1 = 4./3.*((1. - 1.8*dMu*dMu)*sqrt(1. + dLdM*dLdM - dLdM));
1242  k1 = 4./3.*((1. - 1.8*dMu*dMu)*sqrt(1. + dLdM*dLdM) - dLdM);
1243  k2 = -2.*dMu;
1244  } break;
1245 #endif
1246 
1247  default:
1248  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1249  }
1250 
1251  doublereal dd = 1. + dr*(k1*cos(dp) + k2*sin(dp));
1252 
1253  return RRot3*(dd*dUMeanPrev);
1254 };
1255 
1256 /* GlauertRotor - end */
1257 
1258 
1259 /* ManglerRotor - begin */
1260 
1261 ManglerRotor::ManglerRotor(unsigned int uLabel, const DofOwner* pDO)
1262 : Elem(uLabel, flag(0)),
1263 Rotor(uLabel, pDO)
1264 {
1265  NO_OP;
1266 }
1267 
1268 ManglerRotor::ManglerRotor(unsigned int uLabel,
1269  const DofOwner* pDO,
1270  const StructNode* pCraft,
1271  const Mat3x3& rrot,
1272  const StructNode* pRotor,
1273  const StructNode* pGround,
1274  ResForceSet **ppres,
1275  const doublereal& dOR,
1276  const doublereal& dR,
1277  DriveCaller *pdW,
1278  unsigned int iMaxIt,
1279  const doublereal& dTol,
1280  const doublereal& dE,
1281  const doublereal& dCH,
1282  const doublereal& dCFF,
1283  flag fOut)
1284 : Elem(uLabel, flag(0)),
1285 Rotor(uLabel, pDO)
1286 {
1287  Init(pCraft, rrot, pRotor, pGround, ppres, dOR, dR, pdW, iMaxIt, dTol, dE, dCH, dCFF, fOut);
1288 }
1289 
1290 void
1291 ManglerRotor::Init(const StructNode* pCraft,
1292  const Mat3x3& rrot,
1293  const StructNode* pRotor,
1294  const StructNode* pGround,
1295  ResForceSet **ppres,
1296  const doublereal& dOR,
1297  const doublereal& dR,
1298  DriveCaller *pdW,
1299  unsigned int iMaxIt,
1300  const doublereal& dTol,
1301  const doublereal& dE,
1302  const doublereal& dCH,
1303  const doublereal& dCFF,
1304  flag fOut)
1305 {
1306  ASSERT(dOR > 0.);
1307  ASSERT(dR > 0.);
1308  ASSERT(pdW != 0);
1309 
1310  Rotor::Init(pCraft, rrot, pRotor, pGround, ppres, dR, iMaxIt, dTol, dE, fOut);
1311 
1312  dOmegaRef = dOR;
1313  dVTipRef = dOmegaRef*dRadius;
1314  dArea = M_PI*dRadius*dRadius;
1315  Weight.Set(pdW);
1316  dHoverCorrection = dCH;
1317  dForwardFlightCorrection = dCFF;
1318 
1319 #ifdef USE_MPI
1320  if (is_parallel) {
1321  SAFENEWARR(pBlockLenght, int, 18);
1322  SAFENEWARR(pDispl, MPI::Aint, 18);
1323  for (int i = 0; i < 18; i++) {
1324  pBlockLenght[i] = 1;
1325  }
1326  for (int i = 0; i < 3; i++) {
1327  pDispl[i] = MPI::Get_address(&(RRot3.pGetVec()[i]));
1328  }
1329  pDispl[3] = MPI::Get_address(&dUMeanPrev);
1330  pDispl[4] = MPI::Get_address(&dSinAlphad);
1331  pDispl[5] = MPI::Get_address(&dPsi0);
1332  for (int i = 6; i <= 8; i++) {
1333  pDispl[i] = MPI::Get_address(&(Res.Pole().pGetVec()[i-6]));
1334  }
1335  for (int i = 9; i < 18; i++) {
1336  pDispl[i] = MPI::Get_address(&(RRotTranspose.pGetMat()[i-9]));
1337  }
1338  SAFENEWWITHCONSTRUCTOR(pIndVelDataType, MPI::Datatype,
1339  MPI::Datatype(MPI::DOUBLE.Create_hindexed(18, pBlockLenght, pDispl)));
1340  pIndVelDataType->Commit();
1341  }
1342 #endif /* USE_MPI */
1343 }
1344 
1345 
1346 ManglerRotor::~ManglerRotor(void)
1347 {
1348 #ifdef USE_MPI
1349  SAFEDELETEARR(pBlockLenght);
1350  SAFEDELETEARR(pDispl);
1351  SAFEDELETE(pIndVelDataType);
1352 #endif /* USE_MPI */
1353 }
1354 
1355 
1356 /* assemblaggio residuo */
1357 SubVectorHandler&
1358 ManglerRotor::AssRes(SubVectorHandler& WorkVec,
1359  doublereal /* dCoef */ ,
1360  const VectorHandler& /* XCurr */ ,
1361  const VectorHandler& /* XPrimeCurr */ )
1362 {
1363  DEBUGCOUT("Entering ManglerRotor::AssRes()" << std::endl);
1364 
1365 #ifdef USE_MPI
1366  ExchangeLoads(fToBeOutput());
1367  if (!is_parallel || IndVelComm.Get_rank() == 0)
1368 #endif /* USE_MPI */
1369  {
1370  /* Calcola parametri vari */
1371  Rotor::InitParam();
1372 
1373 #ifdef DEBUG
1374  // Prova:
1375 #if 0
1376  Vec3 XTmp(2.,2.,0.);
1377  doublereal dPsiTmp = dGetPsi(XTmp);
1378  doublereal dXTmp = dGetPos(XTmp);
1379  Vec3 IndV = GetInducedVelocity(XTmp);
1380  std::cout
1381  << "X rotore: " << pRotor->GetXCurr() << std::endl
1382  << "V rotore: " << VCraft << std::endl
1383  << "X punto: " << XTmp << std::endl
1384  << "Omega: " << dOmega << std::endl
1385  << "Velocita': " << dVelocity << std::endl
1386  << "Psi0: " << dPsi0 << std::endl
1387  << "Psi punto: " << dPsiTmp << std::endl
1388  << "Raggio: " << dRadius << std::endl
1389  << "r punto: " << dXTmp << std::endl
1390  << "mu: " << dMu << std::endl
1391  << "lambda: " << dLambda << std::endl
1392  << "cos(ad): " << dCosAlphad << std::endl
1393  << "sin(ad): " << dSinAlphad << std::endl
1394  << "UMean: " << dUMean << std::endl
1395  << "iv punto: " << IndV << std::endl;
1396 #endif
1397 #endif /* DEBUG */
1398  }
1399 
1400 #ifdef USE_MPI
1401  ExchangeVelocity();
1402 #endif /* USE_MPI */
1403 
1404  ResetForce();
1405 
1406  /* Non tocca il residuo */
1407  WorkVec.Resize(0);
1408 
1409 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1410  Done();
1411 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
1412 
1413  return WorkVec;
1414 }
1415 
1416 std::ostream&
1417 ManglerRotor::Restart(std::ostream& out) const
1418 {
1419  return Rotor::Restart(out) << "Mangler, " << dRadius << ", ",
1420  Weight.pGetDriveCaller()->Restart(out)
1421  << ", correction, " << dHoverCorrection
1422  << ", " << dForwardFlightCorrection << ';' << std::endl;
1423 }
1424 
1425 /* Somma alla trazione il contributo di forza di un elemento generico */
1426 void
1427 ManglerRotor::AddForce(const Elem *pEl, const StructNode *pNode,
1428  const Vec3& F, const Vec3& M, const Vec3& X)
1429 {
1430 #ifdef USE_MPI
1431  if (ReqV != MPI::REQUEST_NULL) {
1432  while (!ReqV.Test()) {
1433  MYSLEEP(mysleeptime);
1434  }
1435  }
1436 #endif /* USE_MPI */
1437 
1438 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1439  pthread_mutex_lock(&forces_mutex);
1440  Wait();
1441 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
1442 
1443  /* Solo se deve fare l'output calcola anche il momento */
1444  if (bToBeOutput()) {
1445  Res.AddForces(F, M, X);
1446  InducedVelocity::AddForce(pEl, pNode, F, M, X);
1447  } else {
1448  Res.AddForce(F);
1449  }
1450 
1451 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1452  pthread_mutex_unlock(&forces_mutex);
1453 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
1454 }
1455 
1456 
1457 /* Restituisce ad un elemento la velocita' indotta in base alla posizione
1458  * azimuthale */
1459 Vec3
1460 ManglerRotor::GetInducedVelocity(Elem::Type type,
1461  unsigned uLabel, unsigned uPnt, const Vec3& X) const
1462 {
1463  if (dUMeanPrev == 0.) {
1464  return ::Zero3;
1465  }
1466 
1467 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1468  Wait();
1469 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
1470 
1471  doublereal dr, dp;
1472  GetPos(X, dr, dp);
1473 
1474  doublereal dr2 = dr*dr;
1475  doublereal dm2 = 1. - dr2;
1476  doublereal dm = 0.;
1477  if (dm2 > 0.) {
1478  dm = sqrt(dm2);
1479  }
1480  doublereal da = 1. + std::abs(dSinAlphad);
1481  if (da > 1.e-9) {
1482  da = (1. - std::abs(dSinAlphad))/da;
1483  }
1484  if (da > 0.) {
1485  da = sqrt(da);
1486  }
1487 
1488  // c_0
1489  doublereal dd = 15./4.*dm*dr2;
1490 
1491  // c_1
1492  doublereal dc = -15./256.*M_PI*(9.*dr2 - 4.)*dr*da;
1493  dd -= 4.*dc*cos(dp);
1494 
1495  // c_3
1496  dc = 45./256.*M_PI*pow(da*dr, 3);
1497  dd -= 4.*dc*cos(3.*dp);
1498 
1499  // c_[5:2:infty] = 0
1500 
1501  // c_[2:2:infty] (truncated at 10)
1502  for (int i = 2; i <= 10; i += 2) {
1503  dc = pow(-1., i/2 - 1)*15./8.
1504  *((dm + i)/(i*i - 1.)*(3. - 9.*dr2 + i*i) + 3.*dm)/(i*i - 9.)
1505  *pow(((1. - dm)/(1. + dm))*da, i/2.);
1506  dd -= 4.*dc*cos(i*dp);
1507  }
1508 
1509  return RRot3*(dd*dUMeanPrev);
1510 };
1511 
1512 /* ManglerRotor - end */
1513 
1514 
1515 
1516 /*
1517  * According to
1518  * Dale M. Pitt, David A. Peters,
1519  * Theoretical Prediction of Dynamic-Inflow Derivatives,
1520  * Vertica, 5, 1981, pp. 21-34
1521  *
1522  * 8/(3*pi): uncorrected value
1523  * 128/(75*pi): corrected value
1524  */
1525 
1526 #if 0
1527 static const doublereal dM11 = 8./(3.*M_PI);
1528 #else
1529 static const doublereal dM11 = 128./(75.*M_PI);
1530 #endif
1531 
1532 static const doublereal dM22 = -16./(45.*M_PI);
1533 static const doublereal dM33 = -16./(45.*M_PI);
1534 
1535 /* DynamicInflowRotor - begin */
1536 
1537 DynamicInflowRotor::DynamicInflowRotor(unsigned int uLabel, const DofOwner* pDO)
1538 : Elem(uLabel, flag(0)),
1539 Rotor(uLabel, pDO),
1540 dVConst(0), dVSine(0), dVCosine(0),
1541 dL11(0.), dL13(0.), dL22(0.), dL31(0.), dL33(0.)
1542 {
1543  NO_OP;
1544 }
1545 
1546 DynamicInflowRotor::DynamicInflowRotor(unsigned int uLabel,
1547  const DofOwner* pDO,
1548  const StructNode* pCraft,
1549  const Mat3x3& rrot,
1550  const StructNode* pRotor,
1551  const StructNode* pGround,
1552  ResForceSet **ppres,
1553  const doublereal& dOR,
1554  const doublereal& dR,
1555  unsigned int iMaxIt,
1556  const doublereal& dTol,
1557  const doublereal& dE,
1558  const doublereal& dCH,
1559  const doublereal& dCFF,
1560  const doublereal& dVConstTmp,
1561  const doublereal& dVSineTmp,
1562  const doublereal& dVCosineTmp,
1563  flag fOut)
1564 : Elem(uLabel, flag(0)),
1565 Rotor(uLabel, pDO),
1566 dVConst(0), dVSine(0), dVCosine(0),
1567 dL11(0.), dL13(0.), dL22(0.), dL31(0.), dL33(0.)
1568 {
1569  Init(pCraft, rrot, pRotor, pGround, ppres, dOR, dR,
1570  iMaxIt, dTol, dE, dCH, dCFF,
1571  dVConstTmp, dVSineTmp, dVCosineTmp, fOut);
1572 }
1573 
1574 void
1575 DynamicInflowRotor::Init(const StructNode* pCraft,
1576  const Mat3x3& rrot,
1577  const StructNode* pRotor,
1578  const StructNode* pGround,
1579  ResForceSet **ppres,
1580  const doublereal& dOR,
1581  const doublereal& dR,
1582  unsigned int iMaxIt,
1583  const doublereal& dTol,
1584  const doublereal& dE,
1585  const doublereal& dCH,
1586  const doublereal& dCFF,
1587  const doublereal& dVConstTmp,
1588  const doublereal& dVSineTmp,
1589  const doublereal& dVCosineTmp,
1590  flag fOut)
1591 {
1592  ASSERT(dOR > 0.);
1593  ASSERT(dR > 0.);
1594 
1595  Rotor::Init(pCraft, rrot, pRotor, pGround, ppres, dR, iMaxIt, dTol, dE, fOut);
1596 
1597  dVConst = dVConstTmp;
1598  dVSine = dVSineTmp;
1599  dVCosine = dVCosineTmp;
1600 
1601  dOmegaRef = dOR;
1602  dVTipRef = dOmegaRef*dRadius;
1603  dArea = M_PI*dRadius*dRadius;
1604 
1605  dHoverCorrection = dCH;
1606  dForwardFlightCorrection = dCFF;
1607 
1608  /* Significa che valuta la velocita' indotta media al passo corrente */
1609  dWeight = 0.;
1610 
1611 #ifdef USE_MPI
1612  if (is_parallel) {
1613  SAFENEWARR(pBlockLenght, int, 20);
1614  SAFENEWARR(pDispl, MPI::Aint, 20);
1615  for (int i = 0; i < 20; i++) {
1616  pBlockLenght[i] = 1;
1617  }
1618  for (int i = 0; i < 3; i++) {
1619  pDispl[i] = MPI::Get_address(RRot3.pGetVec()+i);
1620  }
1621  pDispl[3] = MPI::Get_address(&dVConst);
1622  pDispl[4] = MPI::Get_address(&dVSine);
1623  pDispl[5] = MPI::Get_address(&dVCosine);
1624  pDispl[6] = MPI::Get_address(&dOmega);
1625  pDispl[7] = MPI::Get_address(&dPsi0);
1626  for (int i = 8; i <= 10; i++) {
1627  pDispl[i] = MPI::Get_address(Res.Pole().pGetVec()+i-8);
1628  }
1629  for (int i = 11; i < 20; i++) {
1630  pDispl[i] = MPI::Get_address(RRotTranspose.pGetMat()+i-11);
1631  }
1632  SAFENEWWITHCONSTRUCTOR(pIndVelDataType, MPI::Datatype,
1633  MPI::Datatype(MPI::DOUBLE.Create_hindexed(20, pBlockLenght, pDispl)));
1634  pIndVelDataType->Commit();
1635  }
1636 #endif /* USE_MPI */
1637 }
1638 
1639 
1640 DynamicInflowRotor::~DynamicInflowRotor(void)
1641 {
1642 #ifdef USE_MPI
1643  SAFEDELETEARR(pBlockLenght);
1644  SAFEDELETEARR(pDispl);
1645  SAFEDELETE(pIndVelDataType);
1646 #endif /* USE_MPI */
1647 }
1648 
1649 
1650 void
1651 DynamicInflowRotor::Output(OutputHandler& OH) const
1652 {
1653  /*
1654  * FIXME: posso usare dei temporanei per il calcolo della trazione
1655  * totale per l'output, cosi' evito il giro dei cast
1656  */
1657  if (bToBeOutput()) {
1658 #ifdef USE_MPI
1659  if (is_parallel && IndVelComm.Get_size() > 1) {
1660  if (IndVelComm.Get_rank() == 0) {
1661  Vec3 TmpF(pTmpVecR), TmpM(pTmpVecR+3);
1662 
1663  OH.Rotors()
1664  << std::setw(8) << GetLabel() /* 1 */
1665  << " " << RRotTranspose*TmpF /* 2-4 */
1666  << " " << RRotTranspose*TmpM /* 5-7 */
1667  << " " << dUMean /* 8 */
1668  << " " << dVelocity /* 9 */
1669  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
1670  << " " << dMu /* 11 */
1671  << " " << dLambda /* 12 */
1672  << " " << dChi /* 13 */
1673  << " " << dPsi0 /* 14 */
1674  << " " << bUMeanRefConverged /* 15 */
1675  << " " << iCurrIter /* 16 */
1676  << " " << dVConst /* 17 */
1677  << " " << dVSine /* 18 */
1678  << " " << dVCosine /* 19 */
1679  << std::endl;
1680 
1681  for (int i = 0; ppRes && ppRes[i]; i++) {
1682  Vec3 TmpF(pTmpVecR+6+6*i);
1683  Vec3 TmpM(pTmpVecR+9+6*i);
1684 
1685  OH.Rotors()
1686  << std::setw(8) << GetLabel()
1687  << ":" << ppRes[i]->GetLabel()
1688  << " " << TmpF
1689  << " " << TmpM
1690  << std::endl;
1691  }
1692  }
1693  } else {
1694  OH.Rotors()
1695  << std::setw(8) << GetLabel() /* 1 */
1696  << " " << RRotTranspose*Res.Force() /* 2-4 */
1697  << " " << RRotTranspose*Res.Moment() /* 5-7 */
1698  << " " << dUMean /* 8 */
1699  << " " << dVelocity /* 9 */
1700  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
1701  << " " << dMu /* 11 */
1702  << " " << dLambda /* 12 */
1703  << " " << dChi /* 13 */
1704  << " " << dPsi0 /* 14 */
1705  << " " << bUMeanRefConverged /* 15 */
1706  << " " << iCurrIter /* 16 */
1707  << " " << dVConst /* 17 */
1708  << " " << dVSine /* 18 */
1709  << " " << dVCosine /* 19 */
1710  << std::endl;
1711 
1712  for (int i = 0; ppRes && ppRes[i]; i++) {
1713  OH.Rotors()
1714  << std::setw(8) << GetLabel()
1715  << ":" << ppRes[i]->GetLabel()
1716  << " " << ppRes[i]->pRes->Force()
1717  << " " << ppRes[i]->pRes->Moment()
1718  << std::endl;
1719  }
1720  }
1721 
1722 #else /* !USE_MPI */
1723  OH.Rotors()
1724  << std::setw(8) << GetLabel() /* 1 */
1725  << " " << RRotTranspose*Res.Force() /* 2-4 */
1726  << " " << RRotTranspose*Res.Moment() /* 5-7 */
1727  << " " << dUMean /* 8 */
1728  << " " << dVelocity /* 9 */
1729  << " " << atan2(dSinAlphad,dCosAlphad) /* 10 */
1730  << " " << dMu /* 11 */
1731  << " " << dLambda /* 12 */
1732  << " " << dChi /* 13 */
1733  << " " << dPsi0 /* 14 */
1734  << " " << bUMeanRefConverged /* 15 */
1735  << " " << iCurrIter /* 16 */
1736  << " " << dVConst /* 17 */
1737  << " " << dVSine /* 18 */
1738  << " " << dVCosine /* 19 */
1739  << std::endl;
1740 
1741  for (int i = 0; ppRes && ppRes[i]; i++) {
1742  OH.Rotors()
1743  << std::setw(8) << GetLabel()
1744  << ":" << ppRes[i]->GetLabel()
1745  << " " << ppRes[i]->pRes->Force()
1746  << " " << ppRes[i]->pRes->Moment()
1747  << std::endl;
1748  }
1749 #endif /* !USE_MPI */
1750  }
1751 }
1752 
1753 /* assemblaggio jacobiano */
1754 VariableSubMatrixHandler&
1755 DynamicInflowRotor::AssJac(VariableSubMatrixHandler& WorkMat,
1756  doublereal dCoef,
1757  const VectorHandler& /* XCurr */ ,
1758  const VectorHandler& /* XPrimeCurr */ )
1759 {
1760  DEBUGCOUT("Entering DynamicInflowRotor::AssJac()" << std::endl);
1761 
1762  WorkMat.SetNullMatrix();
1763 
1764 #ifdef USE_MPI
1765  if (is_parallel && IndVelComm.Get_rank() == 0)
1766 #endif /* USE_MPI */
1767  {
1768  SparseSubMatrixHandler& WM = WorkMat.SetSparse();
1769  integer iFirstIndex = iGetFirstIndex();
1770 
1771  WM.ResizeReset(5, 0);
1772 
1773  WM.PutItem(1, iFirstIndex + 1, iFirstIndex + 1, dM11 + dCoef*dL11);
1774  WM.PutItem(2, iFirstIndex + 3, iFirstIndex + 1, dCoef*dL31);
1775  WM.PutItem(3, iFirstIndex + 2, iFirstIndex + 2, dM22 + dCoef*dL22);
1776  WM.PutItem(4, iFirstIndex + 1, iFirstIndex + 3, dCoef*dL13);
1777  WM.PutItem(5, iFirstIndex + 3, iFirstIndex + 3, dM33 + dCoef*dL33);
1778  }
1779 
1780  return WorkMat;
1781 }
1782 
1783 /* assemblaggio residuo */
1784 SubVectorHandler&
1785 DynamicInflowRotor::AssRes(SubVectorHandler& WorkVec,
1786  doublereal /* dCoef */ ,
1787  const VectorHandler& XCurr,
1788  const VectorHandler& XPrimeCurr)
1789 {
1790  DEBUGCOUT("Entering DynamicInflowRotor::AssRes()" << std::endl);
1791 
1792 #ifdef USE_MPI
1793  ExchangeLoads(flag(1));
1794 
1795  if (!is_parallel || IndVelComm.Get_rank() == 0)
1796 #endif /* USE_MPI */
1797  {
1798  /* Calcola parametri vari */
1799  Rotor::InitParam();
1800 
1801  WorkVec.Resize(3);
1802 
1803  integer iFirstIndex = iGetFirstIndex();
1804 
1805  WorkVec.PutRowIndex(1, iFirstIndex + 1);
1806  WorkVec.PutRowIndex(2, iFirstIndex + 2);
1807  WorkVec.PutRowIndex(3, iFirstIndex + 3);
1808 
1809  dVConst = XCurr(iFirstIndex + 1);
1810  dVSine = XCurr(iFirstIndex + 2);
1811  dVCosine = XCurr(iFirstIndex + 3);
1812 
1813  doublereal dVConstPrime = XPrimeCurr(iFirstIndex + 1);
1814  doublereal dVSinePrime = XPrimeCurr(iFirstIndex + 2);
1815  doublereal dVCosinePrime = XPrimeCurr(iFirstIndex + 3);
1816 
1817  doublereal dCT = 0.;
1818  doublereal dCl = 0.;
1819  doublereal dCm = 0.;
1820 
1821  /*
1822  * Attenzione: moltiplico tutte le equazioni per dOmega
1823  * (ovvero, i coefficienti CT, CL e CM sono divisi
1824  * per dOmega anziche' dOmega^2)
1825  */
1826 
1827  dL11 = 0.;
1828  dL13 = 0.;
1829  dL22 = 0.;
1830  dL31 = 0.;
1831  dL33 = 0.;
1832 
1833  doublereal dDim = dGetAirDensity(GetXCurr())*dArea*dOmega*(dRadius*dRadius);
1834  if (dDim > std::numeric_limits<doublereal>::epsilon()) {
1835  /*
1836  * From Claudio Monteggia:
1837  *
1838  * Ct
1839  * Um = -------------------------------------
1840  * sqrt( lambda^2 / KH^4 + mu^2 / KF^2)
1841  */
1842  doublereal dLambdaTmp
1843  = dLambda/(dHoverCorrection*dHoverCorrection);
1844  doublereal dMuTmp = dMu/dForwardFlightCorrection;
1845 
1846  doublereal dVT
1847  = sqrt(dLambdaTmp*dLambdaTmp + dMuTmp*dMuTmp);
1848  doublereal dVm = 0.;
1849  if (dVT > dVTipTreshold*dVTipRef) {
1850  dVm = (dMuTmp*dMuTmp + dLambdaTmp*(dLambdaTmp + dVConst))/dVT;
1851  }
1852 
1853  /*
1854  * dUMean is just for output;
1855  */
1856  dUMean = dVConst*dOmega*dRadius;
1857 
1858  /* Trazione nel sistema rotore */
1859  doublereal dT = RRot3*Res.Force();
1860 
1861  /* Momento nel sistema rotore-vento */
1862  doublereal dCosP = cos(dPsi0);
1863  doublereal dSinP = sin(dPsi0);
1864  Mat3x3 RTmp( dCosP, -dSinP, 0.,
1865  dSinP, dCosP, 0.,
1866  0., 0., 1.);
1867 
1868  Vec3 M(RTmp*(RRotTranspose*Res.Moment()));
1869 
1870  /* Thrust, roll and pitch coefficients */
1871  dCT = dT/dDim;
1872  dDim *= dRadius;
1873  dCl = - M(1)/dDim;
1874  dCm = - M(2)/dDim;
1875 
1876  if (dVT > std::numeric_limits<doublereal>::epsilon()
1877  && dVm > std::numeric_limits<doublereal>::epsilon())
1878  {
1879 
1880  /* Matrix coefficients */
1881  /* FIXME: divide by 0? */
1882  doublereal dl11 = .5/dVT;
1883  /* FIXME: divide by 0? */
1884  doublereal d = 15./64.*M_PI*tan(dChi/2.);
1885  /* FIXME: divide by 0? */
1886  doublereal dl13 = d/dVm;
1887  /* FIXME: divide by 0? */
1888  doublereal dl31 = d/dVT;
1889 
1890  doublereal dCosChi2 = cos(dChi/2.);
1891  d = 2.*dCosChi2*dCosChi2;
1892  /* FIXME: divide by 0? */
1893  doublereal dl22 = -4./(d*dVm);
1894  /* FIXME: divide by 0? */
1895  doublereal dl33 = -4.*(d - 1)/(d*dVm);
1896 
1897  d = dl11*dl33 - dl31*dl13;
1898  /* FIXME: divide by 0? */
1899  dL11 = dOmega*dl33/d;
1900  dL31 = -dOmega*dl31/d;
1901  dL13 = -dOmega*dl13/d;
1902  dL33 = dOmega*dl11/d;
1903  dL22 = dOmega/dl22;
1904  }
1905  }
1906 
1907 #ifdef DEBUG
1908  /* Prova: */
1909  static int i = -1;
1910  int iv[] = { 0, 1, 0, -1, 0 };
1911  if (++i == 4) {
1912  i = 0;
1913  }
1914  Vec3 XTmp(pRotor->GetXCurr()+pCraft->GetRCurr()*Vec3(dRadius*iv[i],dRadius*iv[i+1],0.));
1915  doublereal dPsiTmp, dXTmp;
1916  GetPos(XTmp, dXTmp, dPsiTmp);
1917  Vec3 IndV = GetInducedVelocity(GetElemType(), GetLabel(), 0, XTmp);
1918  std::cout
1919  << "X rotore: " << pRotor->GetXCurr() << std::endl
1920  << "V rotore: " << VCraft << std::endl
1921  << "X punto: " << XTmp << std::endl
1922  << "Omega: " << dOmega << std::endl
1923  << "Velocita': " << dVelocity << std::endl
1924  << "Psi0: " << dPsi0 << " ("
1925  << dPsi0*180./M_PI << " deg)" << std::endl
1926  << "Psi punto: " << dPsiTmp << " ("
1927  << dPsiTmp*180./M_PI << " deg)" << std::endl
1928  << "Raggio: " << dRadius << std::endl
1929  << "r punto: " << dXTmp << std::endl
1930  << "mu: " << dMu << std::endl
1931  << "lambda: " << dLambda << std::endl
1932  << "cos(ad): " << dCosAlphad << std::endl
1933  << "sin(ad): " << dSinAlphad << std::endl
1934  << "UMean: " << dUMean << std::endl
1935  << "iv punto: " << IndV << std::endl;
1936 #endif /* DEBUG */
1937 
1938  WorkVec.PutCoef(1, dCT - dM11*dVConstPrime
1939  - dL11*dVConst - dL13*dVCosine);
1940  WorkVec.PutCoef(2, dCl - dM22*dVSinePrime
1941  - dL22*dVSine);
1942  WorkVec.PutCoef(3, dCm - dM33*dVCosinePrime
1943  - dL31*dVConst - dL33*dVCosine);
1944 
1945 #ifdef USE_MPI
1946  } else {
1947  WorkVec.Resize(0);
1948 #endif /* USE_MPI */
1949  }
1950 
1951 #ifdef USE_MPI
1952  ExchangeVelocity();
1953 #endif /* USE_MPI */
1954 
1955  /* Ora la trazione non serve piu' */
1956  ResetForce();
1957 
1958 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
1959  Done();
1960 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
1961 
1962  return WorkVec;
1963 }
1964 
1965 /* Relativo ai ...WithDofs */
1966 void
1967 DynamicInflowRotor::SetInitialValue(VectorHandler& /* X */ )
1968 {
1969  NO_OP;
1970 }
1971 
1972 
1973 /* Relativo ai ...WithDofs */
1974 void
1975 DynamicInflowRotor::SetValue(DataManager *pDM,
1976  VectorHandler& X, VectorHandler& XP,
1977  SimulationEntity::Hints *ph)
1978 {
1979  integer iFirstIndex = iGetFirstIndex();
1980 
1981  for (int iCnt = 1; iCnt <= 3; iCnt++) {
1982  XP.PutCoef(iFirstIndex + iCnt, 0.);
1983  }
1984 
1985  X.PutCoef(iFirstIndex + 1, dVConst);
1986  X.PutCoef(iFirstIndex + 2, dVSine);
1987  X.PutCoef(iFirstIndex + 3, dVCosine);
1988 }
1989 
1990 
1991 /* Restart */
1992 std::ostream&
1993 DynamicInflowRotor::Restart(std::ostream& out) const
1994 {
1995  return Rotor::Restart(out) << "dynamic inflow, " << dRadius
1996  << ", correction, " << dHoverCorrection
1997  << ", " << dForwardFlightCorrection << ';' << std::endl;
1998 }
1999 
2000 
2001 /* Somma alla trazione il contributo di forza di un elemento generico */
2002 void
2003 DynamicInflowRotor::AddForce(const Elem *pEl, const StructNode *pNode,
2004  const Vec3& F, const Vec3& M, const Vec3& X)
2005 {
2006  /*
2007  * Gli serve la trazione ed il momento rispetto al rotore,
2008  * che si calcola da se'
2009  */
2010 #ifdef USE_MPI
2011  if (ReqV != MPI::REQUEST_NULL) {
2012  while (!ReqV.Test()) {
2013  MYSLEEP(mysleeptime);
2014  }
2015  }
2016 #endif /* USE_MPI */
2017 
2018 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2019  pthread_mutex_lock(&forces_mutex);
2020  Wait();
2021 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
2022 
2023  Res.AddForces(F, M, X);
2024  if (bToBeOutput()) {
2025  InducedVelocity::AddForce(pEl, pNode, F, M, X);
2026  }
2027 
2028 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2029  pthread_mutex_unlock(&forces_mutex);
2030 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
2031 }
2032 
2033 
2034 /* Restituisce ad un elemento la velocita' indotta in base alla posizione
2035  * azimuthale */
2036 Vec3
2037 DynamicInflowRotor::GetInducedVelocity(Elem::Type type,
2038  unsigned uLabel, unsigned uPnt, const Vec3& X) const
2039 {
2040 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2041  Wait();
2042 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
2043 
2044  doublereal dr, dp;
2045  GetPos(X, dr, dp);
2046 
2047  return RRot3*((dVConst + dr*(dVCosine*cos(dp) + dVSine*sin(dp)))*dRadius*dOmega);
2048 };
2049 
2050 /* DynamicInflowRotor - end */
2051 
2052 
2053 /* PetersHeRotor - begin */
2054 
2055 PetersHeRotor::PetersHeRotor(unsigned int uLabel, const DofOwner* pDO)
2056 : Elem(uLabel, flag(0)),
2057 Rotor(uLabel, pDO),
2058 dVConst(0), dVSine(0), dVCosine(0),
2059 dL11(0.), dL13(0.), dL22(0.), dL31(0.), dL33(0.)
2060 {
2061  NO_OP;
2062 }
2063 
2064 PetersHeRotor::PetersHeRotor(unsigned int uLabel,
2065  const DofOwner* pDO,
2066  const StructNode* pCraft,
2067  const Mat3x3& rrot,
2068  const StructNode* pRotor,
2069  const StructNode* pGround,
2070  ResForceSet **ppres,
2071  const doublereal& dOR,
2072  const doublereal& dR,
2073  unsigned int iMaxIt,
2074  const doublereal& dTol,
2075  const doublereal& dE,
2076  const doublereal& dCH,
2077  const doublereal& dCFF,
2078  const doublereal& dVConstTmp,
2079  const doublereal& dVSineTmp,
2080  const doublereal& dVCosineTmp,
2081  flag fOut)
2082 : Elem(uLabel, flag(0)),
2083 Rotor(uLabel, pDO),
2084 dVConst(0), dVSine(0), dVCosine(0),
2085 dL11(0.), dL13(0.), dL22(0.), dL31(0.), dL33(0.)
2086 {
2087  Init(pCraft, rrot, pRotor, pGround, ppres, dOR, dR,
2088  iMaxIt, dTol, dE, dCH, dCFF,
2089  dVConstTmp, dVSineTmp, dVCosineTmp, fOut);
2090 }
2091 
2092 void
2093 PetersHeRotor::Init(const StructNode* pCraft,
2094  const Mat3x3& rrot,
2095  const StructNode* pRotor,
2096  const StructNode* pGround,
2097  ResForceSet **ppres,
2098  const doublereal& dOR,
2099  const doublereal& dR,
2100  unsigned int iMaxIt,
2101  const doublereal& dTol,
2102  const doublereal& dE,
2103  const doublereal& dCH,
2104  const doublereal& dCFF,
2105  const doublereal& dVConstTmp,
2106  const doublereal& dVSineTmp,
2107  const doublereal& dVCosineTmp,
2108  flag fOut)
2109 {
2110  ASSERT(dOR > 0.);
2111  ASSERT(dR > 0.);
2112 
2113  Rotor::Init(pCraft, rrot, pRotor, pGround, ppres, dR, iMaxIt, dTol, dE, fOut);
2114 
2115  dVConst = dVConstTmp;
2116  dVSine = dVSineTmp;
2117  dVCosine = dVCosineTmp;
2118 
2119  dOmegaRef = dOR;
2120  dVTipRef = dOmegaRef*dRadius;
2121  dArea = M_PI*dRadius*dRadius;
2122 
2123  dHoverCorrection = dCH;
2124  dForwardFlightCorrection = dCFF;
2125 
2126  /* Significa che valuta la velocita' indotta media al passo corrente */
2127  dWeight = 0.;
2128 
2129 #ifdef USE_MPI
2130  if (is_parallel) {
2131  SAFENEWARR(pBlockLenght, int, 20);
2132  SAFENEWARR(pDispl, MPI::Aint, 20);
2133  for (int i = 0; i < 20; i++) {
2134  pBlockLenght[i] = 1;
2135  }
2136  for (int i = 0; i < 3; i++) {
2137  pDispl[i] = MPI::Get_address(RRot3.pGetVec()+i);
2138  }
2139  pDispl[3] = MPI::Get_address(&dVConst);
2140  pDispl[4] = MPI::Get_address(&dVSine);
2141  pDispl[5] = MPI::Get_address(&dVCosine);
2142  pDispl[6] = MPI::Get_address(&dOmega);
2143  pDispl[7] = MPI::Get_address(&dPsi0);
2144  for (int i = 8; i <= 10; i++) {
2145  pDispl[i] = MPI::Get_address(Res.Pole().pGetVec()+i-8);
2146  }
2147  for (int i = 11; i < 20; i++) {
2148  pDispl[i] = MPI::Get_address(RRotTranspose.pGetMat()+i-11);
2149  }
2150  SAFENEWWITHCONSTRUCTOR(pIndVelDataType, MPI::Datatype,
2151  MPI::Datatype(MPI::DOUBLE.Create_hindexed(20, pBlockLenght, pDispl)));
2152  pIndVelDataType->Commit();
2153  }
2154 #endif /* USE_MPI */
2155 }
2156 
2157 
2158 PetersHeRotor::~PetersHeRotor(void)
2159 {
2160 #ifdef USE_MPI
2161  SAFEDELETEARR(pBlockLenght);
2162  SAFEDELETEARR(pDispl);
2163  SAFEDELETE(pIndVelDataType);
2164 #endif /* USE_MPI */
2165 }
2166 
2167 
2168 void
2169 PetersHeRotor::Output(OutputHandler& OH) const
2170 {
2171  /*
2172  * FIXME: posso usare dei temporanei per il calcolo della trazione
2173  * totale per l'output, cosi' evito il giro dei cast
2174  */
2175  if (bToBeOutput()) {
2176 #ifdef USE_MPI
2177  if (is_parallel && IndVelComm.Get_size() > 1) {
2178  if (IndVelComm.Get_rank() == 0) {
2179  Vec3 TmpF(pTmpVecR), TmpM(pTmpVecR+3);
2180 
2181  OH.Rotors()
2182  << std::setw(8) << GetLabel() /* 1 */
2183  << " " << RRotTranspose*TmpF /* 2-4 */
2184  << " " << RRotTranspose*TmpM /* 5-7 */
2185  << " " << dUMean /* 8 */
2186  << " " << dVelocity /* 9 */
2187  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
2188  << " " << dMu /* 11 */
2189  << " " << dLambda /* 12 */
2190  << " " << dChi /* 13 */
2191  << " " << dPsi0 /* 14 */
2192  << " " << bUMeanRefConverged /* 15 */
2193  << " " << iCurrIter /* 16 */
2194  << " " << dVConst /* 17 */
2195  << " " << dVSine /* 18 */
2196  << " " << dVCosine /* 19 */
2197  << std::endl;
2198 
2199  for (int i = 0; ppRes && ppRes[i]; i++) {
2200  Vec3 TmpF(pTmpVecR+6+6*i);
2201  Vec3 TmpM(pTmpVecR+9+6*i);
2202 
2203  OH.Rotors()
2204  << std::setw(8) << GetLabel()
2205  << ":" << ppRes[i]->GetLabel()
2206  << " " << TmpF
2207  << " " << TmpM
2208  << std::endl;
2209  }
2210  }
2211  } else {
2212  OH.Rotors()
2213  << std::setw(8) << GetLabel() /* 1 */
2214  << " " << RRotTranspose*Res.Force() /* 2-4 */
2215  << " " << RRotTranspose*Res.Moment() /* 5-7 */
2216  << " " << dUMean /* 8 */
2217  << " " << dVelocity /* 9 */
2218  << " " << atan2(dSinAlphad, dCosAlphad) /* 10 */
2219  << " " << dMu /* 11 */
2220  << " " << dLambda /* 12 */
2221  << " " << dChi /* 13 */
2222  << " " << dPsi0 /* 14 */
2223  << " " << bUMeanRefConverged /* 15 */
2224  << " " << iCurrIter /* 16 */
2225  << " " << dVConst /* 17 */
2226  << " " << dVSine /* 18 */
2227  << " " << dVCosine /* 19 */
2228  << std::endl;
2229 
2230  for (int i = 0; ppRes && ppRes[i]; i++) {
2231  OH.Rotors()
2232  << std::setw(8) << GetLabel()
2233  << ":" << ppRes[i]->GetLabel()
2234  << " " << ppRes[i]->pRes->Force()
2235  << " " << ppRes[i]->pRes->Moment()
2236  << std::endl;
2237  }
2238  }
2239 
2240 #else /* !USE_MPI */
2241  OH.Rotors()
2242  << std::setw(8) << GetLabel() /* 1 */
2243  << " " << RRotTranspose*Res.Force() /* 2-4 */
2244  << " " << RRotTranspose*Res.Moment() /* 5-7 */
2245  << " " << dUMean /* 8 */
2246  << " " << dVelocity /* 9 */
2247  << " " << atan2(dSinAlphad,dCosAlphad) /* 10 */
2248  << " " << dMu /* 11 */
2249  << " " << dLambda /* 12 */
2250  << " " << dChi /* 13 */
2251  << " " << dPsi0 /* 14 */
2252  << " " << bUMeanRefConverged /* 15 */
2253  << " " << iCurrIter /* 16 */
2254  << " " << dVConst /* 17 */
2255  << " " << dVSine /* 18 */
2256  << " " << dVCosine /* 19 */
2257  << std::endl;
2258 
2259  for (int i = 0; ppRes && ppRes[i]; i++) {
2260  OH.Rotors()
2261  << std::setw(8) << GetLabel()
2262  << ":" << ppRes[i]->GetLabel()
2263  << " " << ppRes[i]->pRes->Force()
2264  << " " << ppRes[i]->pRes->Moment()
2265  << std::endl;
2266  }
2267 #endif /* !USE_MPI */
2268  }
2269 }
2270 
2271 /* assemblaggio jacobiano */
2272 VariableSubMatrixHandler&
2273 PetersHeRotor::AssJac(VariableSubMatrixHandler& WorkMat,
2274  doublereal dCoef,
2275  const VectorHandler& /* XCurr */ ,
2276  const VectorHandler& /* XPrimeCurr */ )
2277 {
2278  DEBUGCOUT("Entering PetersHeRotor::AssJac()" << std::endl);
2279 
2280  WorkMat.SetNullMatrix();
2281 
2282 #ifdef USE_MPI
2283  if (is_parallel && IndVelComm.Get_rank() == 0)
2284 #endif /* USE_MPI */
2285  {
2286  SparseSubMatrixHandler& WM = WorkMat.SetSparse();
2287  integer iFirstIndex = iGetFirstIndex();
2288 
2289  WM.ResizeReset(5, 0);
2290 
2291  WM.PutItem(1, iFirstIndex + 1, iFirstIndex + 1, dM11 + dCoef*dL11);
2292  WM.PutItem(2, iFirstIndex + 3, iFirstIndex + 1, dCoef*dL31);
2293  WM.PutItem(3, iFirstIndex + 2, iFirstIndex + 2, dM22 + dCoef*dL22);
2294  WM.PutItem(4, iFirstIndex + 1, iFirstIndex + 3, dCoef*dL13);
2295  WM.PutItem(5, iFirstIndex + 3, iFirstIndex + 3, dM33 + dCoef*dL33);
2296  }
2297 
2298  return WorkMat;
2299 }
2300 
2301 /* assemblaggio residuo */
2302 SubVectorHandler&
2303 PetersHeRotor::AssRes(SubVectorHandler& WorkVec,
2304  doublereal /* dCoef */ ,
2305  const VectorHandler& XCurr,
2306  const VectorHandler& XPrimeCurr)
2307 {
2308  DEBUGCOUT("Entering PetersHeRotor::AssRes()" << std::endl);
2309 
2310 #ifdef USE_MPI
2311  ExchangeLoads(flag(1));
2312 
2313  if (!is_parallel || IndVelComm.Get_rank() == 0)
2314 #endif /* USE_MPI */
2315  {
2316  /* Calcola parametri vari */
2317  Rotor::InitParam();
2318 
2319  WorkVec.Resize(3);
2320 
2321  integer iFirstIndex = iGetFirstIndex();
2322 
2323  WorkVec.PutRowIndex(1, iFirstIndex + 1);
2324  WorkVec.PutRowIndex(2, iFirstIndex + 2);
2325  WorkVec.PutRowIndex(3, iFirstIndex + 3);
2326 
2327  dVConst = XCurr(iFirstIndex + 1);
2328  dVSine = XCurr(iFirstIndex + 2);
2329  dVCosine = XCurr(iFirstIndex + 3);
2330 
2331  doublereal dVConstPrime = XPrimeCurr(iFirstIndex + 1);
2332  doublereal dVSinePrime = XPrimeCurr(iFirstIndex + 2);
2333  doublereal dVCosinePrime = XPrimeCurr(iFirstIndex + 3);
2334 
2335  doublereal dCT = 0.;
2336  doublereal dCl = 0.;
2337  doublereal dCm = 0.;
2338 
2339  /*
2340  * Attenzione: moltiplico tutte le equazioni per dOmega
2341  * (ovvero, i coefficienti CT, CL e CM sono divisi
2342  * per dOmega anziche' dOmega^2)
2343  */
2344 
2345  dL11 = 0.;
2346  dL13 = 0.;
2347  dL22 = 0.;
2348  dL31 = 0.;
2349  dL33 = 0.;
2350 
2351  doublereal dDim = dGetAirDensity(GetXCurr())*dArea*dOmega*(dRadius*dRadius);
2352  if (dDim > std::numeric_limits<doublereal>::epsilon()) {
2353  /*
2354  * From Claudio Monteggia:
2355  *
2356  * Ct
2357  * Um = -------------------------------------
2358  * sqrt( lambda^2 / KH^4 + mu^2 / KF^2)
2359  */
2360  doublereal dLambdaTmp
2361  = dLambda/(dHoverCorrection*dHoverCorrection);
2362  doublereal dMuTmp = dMu/dForwardFlightCorrection;
2363 
2364  doublereal dVT
2365  = sqrt(dLambdaTmp*dLambdaTmp + dMuTmp*dMuTmp);
2366  doublereal dVm = 0.;
2367  if (dVT > dVTipTreshold*dVTipRef) {
2368  dVm = (dMuTmp*dMuTmp + dLambdaTmp*(dLambdaTmp + dVConst))/dVT;
2369  }
2370 
2371  /*
2372  * dUMean is just for output;
2373  */
2374  dUMean = dVConst*dOmega*dRadius;
2375 
2376  /* Trazione nel sistema rotore */
2377  doublereal dT = RRot3*Res.Force();
2378 
2379  /* Momento nel sistema rotore-vento */
2380  doublereal dCosP = cos(dPsi0);
2381  doublereal dSinP = sin(dPsi0);
2382  Mat3x3 RTmp( dCosP, -dSinP, 0.,
2383  dSinP, dCosP, 0.,
2384  0., 0., 1.);
2385 
2386  Vec3 M(RTmp*(RRotTranspose*Res.Moment()));
2387 
2388  /* Thrust, roll and pitch coefficients */
2389  dCT = dT/dDim;
2390  dDim *= dRadius;
2391  dCl = - M(1)/dDim;
2392  dCm = - M(2)/dDim;
2393 
2394  if (dVT > std::numeric_limits<doublereal>::epsilon()
2395  && dVm > std::numeric_limits<doublereal>::epsilon())
2396  {
2397 
2398  /* Matrix coefficients */
2399  /* FIXME: divide by 0? */
2400  doublereal dl11 = .5/dVT;
2401  /* FIXME: divide by 0? */
2402  doublereal d = 15./64.*M_PI*tan(dChi/2.);
2403  /* FIXME: divide by 0? */
2404  doublereal dl13 = d/dVm;
2405  /* FIXME: divide by 0? */
2406  doublereal dl31 = d/dVT;
2407 
2408  doublereal dCosChi2 = cos(dChi/2.);
2409  d = 2.*dCosChi2*dCosChi2;
2410  /* FIXME: divide by 0? */
2411  doublereal dl22 = -4./(d*dVm);
2412  /* FIXME: divide by 0? */
2413  doublereal dl33 = -4.*(d - 1)/(d*dVm);
2414 
2415  d = dl11*dl33 - dl31*dl13;
2416  /* FIXME: divide by 0? */
2417  dL11 = dOmega*dl33/d;
2418  dL31 = -dOmega*dl31/d;
2419  dL13 = -dOmega*dl13/d;
2420  dL33 = dOmega*dl11/d;
2421  dL22 = dOmega/dl22;
2422  }
2423  }
2424 
2425 #ifdef DEBUG
2426  /* Prova: */
2427  static int i = -1;
2428  int iv[] = { 0, 1, 0, -1, 0 };
2429  if (++i == 4) {
2430  i = 0;
2431  }
2432  Vec3 XTmp(pRotor->GetXCurr()+pCraft->GetRCurr()*Vec3(dRadius*iv[i],dRadius*iv[i+1],0.));
2433  doublereal dPsiTmp, dXTmp;
2434  GetPos(XTmp, dXTmp, dPsiTmp);
2435  Vec3 IndV = GetInducedVelocity(GetElemType(), GetLabel(), 0, XTmp);
2436  std::cout
2437  << "X rotore: " << pRotor->GetXCurr() << std::endl
2438  << "V rotore: " << VCraft << std::endl
2439  << "X punto: " << XTmp << std::endl
2440  << "Omega: " << dOmega << std::endl
2441  << "Velocita': " << dVelocity << std::endl
2442  << "Psi0: " << dPsi0 << " ("
2443  << dPsi0*180./M_PI << " deg)" << std::endl
2444  << "Psi punto: " << dPsiTmp << " ("
2445  << dPsiTmp*180./M_PI << " deg)" << std::endl
2446  << "Raggio: " << dRadius << std::endl
2447  << "r punto: " << dXTmp << std::endl
2448  << "mu: " << dMu << std::endl
2449  << "lambda: " << dLambda << std::endl
2450  << "cos(ad): " << dCosAlphad << std::endl
2451  << "sin(ad): " << dSinAlphad << std::endl
2452  << "UMean: " << dUMean << std::endl
2453  << "iv punto: " << IndV << std::endl;
2454 #endif /* DEBUG */
2455 
2456  WorkVec.PutCoef(1, dCT - dM11*dVConstPrime
2457  - dL11*dVConst - dL13*dVCosine);
2458  WorkVec.PutCoef(2, dCl - dM22*dVSinePrime
2459  - dL22*dVSine);
2460  WorkVec.PutCoef(3, dCm - dM33*dVCosinePrime
2461  - dL31*dVConst - dL33*dVCosine);
2462 
2463 #ifdef USE_MPI
2464  } else {
2465  WorkVec.Resize(0);
2466 #endif /* USE_MPI */
2467  }
2468 
2469 #ifdef USE_MPI
2470  ExchangeVelocity();
2471 #endif /* USE_MPI */
2472 
2473  /* Ora la trazione non serve piu' */
2474  ResetForce();
2475 
2476 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2477  Done();
2478 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
2479 
2480  return WorkVec;
2481 }
2482 
2483 /* Relativo ai ...WithDofs */
2484 void
2485 PetersHeRotor::SetInitialValue(VectorHandler& /* X */ )
2486 {
2487  NO_OP;
2488 }
2489 
2490 
2491 /* Relativo ai ...WithDofs */
2492 void
2493 PetersHeRotor::SetValue(DataManager *pDM,
2494  VectorHandler& X, VectorHandler& XP,
2495  SimulationEntity::Hints *ph)
2496 {
2497  integer iFirstIndex = iGetFirstIndex();
2498 
2499  for (int iCnt = 1; iCnt <= 3; iCnt++) {
2500  XP.PutCoef(iFirstIndex + iCnt, 0.);
2501  }
2502 
2503  X.PutCoef(iFirstIndex + 1, dVConst);
2504  X.PutCoef(iFirstIndex + 2, dVSine);
2505  X.PutCoef(iFirstIndex + 3, dVCosine);
2506 }
2507 
2508 
2509 /* Restart */
2510 std::ostream&
2511 PetersHeRotor::Restart(std::ostream& out) const
2512 {
2513  return Rotor::Restart(out) << "dynamic inflow, " << dRadius
2514  << ", correction, " << dHoverCorrection
2515  << ", " << dForwardFlightCorrection << ';' << std::endl;
2516 }
2517 
2518 
2519 /* Somma alla trazione il contributo di forza di un elemento generico */
2520 void
2521 PetersHeRotor::AddForce(const Elem *pEl, const StructNode *pNode,
2522  const Vec3& F, const Vec3& M, const Vec3& X)
2523 {
2524  /*
2525  * Gli serve la trazione ed il momento rispetto al rotore,
2526  * che si calcola da se'
2527  */
2528 #ifdef USE_MPI
2529  if (ReqV != MPI::REQUEST_NULL) {
2530  while (!ReqV.Test()) {
2531  MYSLEEP(mysleeptime);
2532  }
2533  }
2534 #endif /* USE_MPI */
2535 
2536 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2537  pthread_mutex_lock(&forces_mutex);
2538  Wait();
2539 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
2540 
2541  Res.AddForces(F, M, X);
2542  if (bToBeOutput()) {
2543  InducedVelocity::AddForce(pEl, pNode, F, M, X);
2544  }
2545 
2546 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2547  pthread_mutex_unlock(&forces_mutex);
2548 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
2549 }
2550 
2551 
2552 /* Restituisce ad un elemento la velocita' indotta in base alla posizione
2553  * azimuthale */
2554 Vec3
2555 PetersHeRotor::GetInducedVelocity(Elem::Type type,
2556  unsigned uLabel, unsigned uPnt, const Vec3& X) const
2557 {
2558 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
2559  Wait();
2560 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
2561 
2562  doublereal dr, dp;
2563  GetPos(X, dr, dp);
2564 
2565  return RRot3*((dVConst + dr*(dVCosine*cos(dp) + dVSine*sin(dp)))*dRadius*dOmega);
2566 };
2567 
2568 /* PetersHeRotor - end */
2569 
2570 
2571 /* Legge un rotore */
2572 Elem *
2573 ReadRotor(DataManager* pDM,
2574  MBDynParser& HP,
2575  const DofOwner* pDO,
2576  unsigned int uLabel)
2577 {
2578  DEBUGCOUT("Entering ReadRotor()" << std::endl);
2579 
2580  /* demote to pedantic; syntax changed a long ago... */
2581  pedantic_cout("WARNING: the syntax changed; use a comma ',' "
2582  "instead of a colon ':' after the keyword "
2583  "\"induced velocity\"" << std::endl);
2584 
2585  const char* sKeyWords[] = {
2586  "induced" "velocity",
2587  "no",
2588  "uniform",
2589  "uniform" "sectional",
2590  "glauert",
2591  "mangler",
2592  "dynamic" "inflow",
2593 
2594  NULL
2595  };
2596 
2597  /* enum delle parole chiave */
2598  enum KeyWords {
2599  UNKNOWN = -1,
2600  INDUCEDVELOCITY = 0,
2601  NO,
2602  UNIFORM,
2603  UNIFORM_SECTIONAL,
2604  GLAUERT,
2605  MANGLER,
2606  DYNAMICINFLOW,
2607  PETERS_HE,
2608 
2609  LASTKEYWORD
2610  };
2611 
2612  /* tabella delle parole chiave */
2613  KeyTable K(HP, sKeyWords);
2614 
2615  /* aircraft node */
2616  const StructNode* pCraft = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
2617 
2618  /* rotor orientation with respect to aircraft */
2619  Mat3x3 rrot(Eye3);
2620  if (HP.IsKeyWord("orientation")) {
2621  ReferenceFrame RF(pCraft);
2622  rrot = HP.GetRotRel(RF);
2623 
2624  } else if (HP.IsKeyWord("hinge")) {
2625  silent_cerr("InducedVelocity(" << uLabel << "): deprecated keyword \"hinge\"; use \"orientation\" instead at line " << HP.GetLineData() << std::endl);
2626 
2627  ReferenceFrame RF(pCraft);
2628  rrot = HP.GetRotRel(RF);
2629  }
2630 
2631  /* rotor node */
2632  const StructNode* pRotor = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
2633 
2634  KeyWords InducedType = NO;
2635  if (HP.IsArg() && HP.IsKeyWord("induced" "velocity")) {
2636  InducedType = KeyWords(HP.GetWord());
2637  }
2638 
2639  Elem* pEl = 0;
2640  ResForceSet **ppres = 0;
2641 
2642  switch (InducedType) {
2643  case NO: {
2644  DEBUGCOUT("No induced velocity is considered" << std::endl);
2645 
2646  doublereal dR = 0.;
2647  if (HP.IsKeyWord("radius")) {
2648  dR = HP.GetReal();
2649  if (dR <= 0) {
2650  silent_cerr("Rotor(" << uLabel << "): "
2651  "invalid radius " << dR
2652  << " at line " << HP.GetLineData()
2653  << std::endl);
2654  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2655  }
2656  }
2657 
2658  ppres = ReadResSets(pDM, HP);
2659 
2660  flag fOut = pDM->fReadOutput(HP, Elem::INDUCEDVELOCITY);
2661 
2662  SAFENEWWITHCONSTRUCTOR(pEl,
2663  NoRotor,
2664  NoRotor(uLabel, pDO, pCraft, rrot, pRotor,
2665  ppres, dR, fOut));
2666  } break;
2667 
2668  case UNIFORM:
2669  case UNIFORM_SECTIONAL:
2670  case GLAUERT:
2671  case MANGLER:
2672  case DYNAMICINFLOW: {
2673  GlauertRotor::Type gtype(GlauertRotor::GLAUERT);
2674  if (InducedType == GLAUERT && HP.IsKeyWord("type")) {
2675  if (HP.IsKeyWord("glauert")) {
2676  gtype = GlauertRotor::GLAUERT;
2677 
2678  } else if (HP.IsKeyWord("coleman")) {
2679  gtype = GlauertRotor::COLEMAN_ET_AL;
2680 
2681  } else if (HP.IsKeyWord("drees")) {
2682  gtype = GlauertRotor::DREES_1;
2683 
2684  } else if (HP.IsKeyWord("payne")) {
2685  gtype = GlauertRotor::PAYNE;
2686 
2687  } else if (HP.IsKeyWord("white" "and" "blake")) {
2688  gtype = GlauertRotor::WHITE_AND_BLAKE;
2689 
2690  } else if (HP.IsKeyWord("pitt" "and" "peters")) {
2691  gtype = GlauertRotor::PITT_AND_PETERS;
2692 
2693  } else if (HP.IsKeyWord("howlett")) {
2694  gtype = GlauertRotor::HOWLETT;
2695 
2696  } else if (HP.IsKeyWord("drees" "2")) {
2697  silent_cerr("warning, \"drees 2\" deprecated at line " << HP.GetLineData() << "; use \"drees\" instead" << std::endl);
2698  gtype = GlauertRotor::DREES_2;
2699 
2700  } else {
2701  silent_cerr("Rotor(" << uLabel << "): "
2702  "unknown variant of Glauert's induced velocity at line "
2703  << HP.GetLineData() << std::endl);
2704  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2705  }
2706  }
2707 
2708  doublereal dOR = HP.GetReal();
2709  DEBUGCOUT("Reference rotation speed: " << dOR << std::endl);
2710  if (dOR <= 0.) {
2711  silent_cerr("Rotor(" << uLabel << "): "
2712  "invalid reference speed " << dOR
2713  << " at line " << HP.GetLineData()
2714  << std::endl);
2715  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2716  }
2717 
2718  doublereal dR = HP.GetReal();
2719  DEBUGCOUT("Radius: " << dR << std::endl);
2720  if (dR <= 0.) {
2721  silent_cerr("Rotor(" << uLabel << "): "
2722  "invalid radius " << dR
2723  << " at line " << HP.GetLineData()
2724  << std::endl);
2725  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
2726  }
2727 
2728  // optional parameters
2729  bool bGotInitialValues(false);
2730  doublereal dVConst = 0.;
2731  doublereal dVSine = 0.;
2732  doublereal dVCosine = 0.;
2733  DriveCaller *pdW = 0;
2734  const StructNode *pGround = 0;
2735  unsigned iMaxIter = unsigned(-1);
2736  doublereal dTolerance = std::numeric_limits<double>::max();
2737  doublereal dEta = -1.;
2738  doublereal dCH = -1.;
2739  doublereal dCFF = -1.;
2740 
2741  while (HP.IsArg()) {
2742  if (HP.IsKeyWord("ground")) {
2743  /*
2744  * ground node for ground effect modeling
2745  */
2746  if (pGround != 0) {
2747  silent_cerr("Rotor(" << uLabel << "): "
2748  "providing another \"ground\" node "
2749  "at line " << HP.GetLineData()
2750  << std::endl);
2751  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2752  }
2753 
2754  /* ground node */
2755  pGround = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
2756 
2757  } else if (HP.IsKeyWord("initial" "value")) {
2758  if (InducedType != DYNAMICINFLOW) {
2759  silent_cerr("Rotor(" << uLabel << "): "
2760  "invalid parameter \"initial value\" "
2761  "at line " << HP.GetLineData()
2762  << std::endl);
2763  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2764  }
2765 
2766  if (bGotInitialValues) {
2767  silent_cerr("Rotor(" << uLabel << "): "
2768  "providing \"initial value\" another time "
2769  "at line " << HP.GetLineData()
2770  << std::endl);
2771  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2772  }
2773 
2774  dVConst = HP.GetReal();
2775  dVSine = HP.GetReal();
2776  dVCosine = HP.GetReal();
2777 
2778  bGotInitialValues = true;
2779 
2780  } else if (HP.IsKeyWord("weight") || HP.IsKeyWord("delay")) {
2781  if (InducedType == DYNAMICINFLOW) {
2782  silent_cerr("Rotor(" << uLabel << "): "
2783  "invalid parameter \"delay\" "
2784  "at line " << HP.GetLineData()
2785  << std::endl);
2786  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2787  }
2788 
2789  if (pdW != 0) {
2790  silent_cerr("Rotor(" << uLabel << "): "
2791  "providing another \"delay\" driver "
2792  "at line " << HP.GetLineData()
2793  << std::endl);
2794  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2795  }
2796 
2797  /*
2798  * Legge il coefficiente di peso della velocita'
2799  * indotta ("weight" e' deprecato, si preferisce
2800  * "delay")
2801  *
2802  * nota:
2803  *
2804  * U = U_n * ( 1 - dW ) + U_n-1 * dW
2805  *
2806  * quindi dW rappresenta il peso che si da'
2807  * al valore al passo precedente; in questo modo
2808  * si introduce un ritardo euristico (attenzione:
2809  * il ritardo vero dipende dal passo temporale)
2810  * che aiuta ad evitare problemi di convergenza.
2811  * Se si desidera un ritardo "fisico", conviene
2812  * provare il "Dynamic Inflow".
2813  */
2814  pdW = HP.GetDriveCaller();
2815 
2816  } else if (HP.IsKeyWord("max" "iterations")) {
2817  if (iMaxIter != unsigned(-1)) {
2818  silent_cerr("Rotor(" << uLabel << "): "
2819  "providing another \"max iterations\" value "
2820  "at line " << HP.GetLineData()
2821  << std::endl);
2822  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2823  }
2824 
2825  /* max iterations when computing reference inflow velocity;
2826  * after iMaxIter iterations, the current value is accepted
2827  * regardless of convergence; thus, 1 reproduces original
2828  * behavior */
2829  int i = HP.GetInt();
2830  if (i <= 0) {
2831  silent_cerr("illegal max iterations "
2832  << i << " for Rotor(" << uLabel << ")");
2833  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2834  }
2835  iMaxIter = i;
2836 
2837  } else if (HP.IsKeyWord("tolerance")) {
2838  if (dTolerance != std::numeric_limits<double>::max()) {
2839  silent_cerr("Rotor(" << uLabel << "): "
2840  "providing another \"tolerance\" value "
2841  "at line " << HP.GetLineData()
2842  << std::endl);
2843  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2844  }
2845 
2846  /* tolerance when computing reference inflow velocity;
2847  * when the difference in inflow velocity between two
2848  * iterations is less than tolerance in module, the
2849  * cycle breaks */
2850 
2851  dTolerance = HP.GetReal();
2852  if (dTolerance <= 0.) {
2853  silent_cerr("illegal tolerance "
2854  << dTolerance << " for Rotor(" << uLabel << ")");
2855  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2856  }
2857 
2858  } else if (HP.IsKeyWord("eta")) {
2859  if (dEta != -1.) {
2860  silent_cerr("Rotor(" << uLabel << "): "
2861  "providing another \"eta\" relaxation factor "
2862  "at line " << HP.GetLineData()
2863  << std::endl);
2864  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2865  }
2866 
2867  /* increment factor when computing reference inflow velocity;
2868  * only a fraction dEta of the difference between two iterations
2869  * is applied */
2870 
2871  dEta = HP.GetReal();
2872  if (dEta <= 0.) {
2873  silent_cerr("illegal eta "
2874  << dEta << " for Rotor(" << uLabel << ")");
2875  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2876  }
2877 
2878  } else if (HP.IsKeyWord("correction")) {
2879  if (dCH != -1.) {
2880  silent_cerr("Rotor(" << uLabel << "): "
2881  "providing another \"correction\" factor "
2882  "at line " << HP.GetLineData()
2883  << std::endl);
2884  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2885  }
2886 
2887  /* Legge la correzione della velocita' indotta */
2888  dCH = HP.GetReal();
2889  DEBUGCOUT("Hover correction: " << dCH << std::endl);
2890  if (dCH <= 0.) {
2891  silent_cerr("Rotor(" << uLabel << "): "
2892  "illegal null or negative hover inflow correction "
2893  "at line " << HP.GetLineData()
2894  << std::endl);
2895  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2896  }
2897 
2898  dCFF = HP.GetReal();
2899  DEBUGCOUT("Forward-flight correction: " << dCFF << std::endl);
2900  if (dCFF <= 0.) {
2901  silent_cerr("Rotor(" << uLabel << "): "
2902  "illegal null or negative forward-flight inflow correction "
2903  "at line " << HP.GetLineData()
2904  << std::endl);
2905  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2906  }
2907 
2908  } else {
2909  break;
2910  }
2911  }
2912 
2913  ppres = ReadResSets(pDM, HP);
2914 
2915  flag fOut = pDM->fReadOutput(HP, Elem::INDUCEDVELOCITY);
2916 
2917  // check consistency and initialize defaults
2918  if (InducedType != DYNAMICINFLOW && pdW == 0) {
2919  SAFENEW(pdW, NullDriveCaller);
2920  }
2921 
2922  if (iMaxIter == unsigned(-1)) {
2923  iMaxIter = 1;
2924 
2925  } else {
2926  if (dTolerance == std::numeric_limits<double>::max()) {
2927  silent_cerr("Rotor(" << uLabel << "): "
2928  "warning, \"max iterations\" is meaningless with default tolerance"
2929  << std::endl);
2930  }
2931  }
2932 
2933  if (dEta == -1.) {
2934  dEta = 1.;
2935  }
2936 
2937  if (dCH == -1.) {
2938  dCH = 1.;
2939  dCFF = 1.;
2940  }
2941 
2942  // create element
2943  switch (InducedType) {
2944  case UNIFORM:
2945  DEBUGCOUT("Uniform induced velocity" << std::endl);
2946  SAFENEWWITHCONSTRUCTOR(pEl,
2947  UniformRotor,
2948  UniformRotor(uLabel, pDO, pCraft, rrot,
2949  pRotor, pGround,
2950  ppres, dOR, dR, pdW,
2951  iMaxIter, dTolerance, dEta,
2952  dCH, dCFF,
2953  fOut));
2954  break;
2955 
2956  case UNIFORM_SECTIONAL:
2957  DEBUGCOUT("Uniform induced velocity" << std::endl);
2958  SAFENEWWITHCONSTRUCTOR(pEl,
2959  UniformRotor2,
2960  UniformRotor2(uLabel, pDO, pCraft, rrot,
2961  pRotor, pGround,
2962  ppres, dOR, dR, pdW,
2963  iMaxIter, dTolerance, dEta,
2964  dCH, dCFF,
2965  fOut));
2966  break;
2967 
2968  case GLAUERT:
2969  DEBUGCOUT("Glauert induced velocity" << std::endl);
2970  SAFENEWWITHCONSTRUCTOR(pEl,
2971  GlauertRotor,
2972  GlauertRotor(uLabel, pDO, pCraft, rrot,
2973  pRotor, pGround,
2974  ppres, dOR, dR, pdW,
2975  iMaxIter, dTolerance, dEta,
2976  dCH, dCFF, gtype,
2977  fOut));
2978  break;
2979 
2980  case MANGLER:
2981  DEBUGCOUT("Mangler induced velocity" << std::endl);
2982 
2983  SAFENEWWITHCONSTRUCTOR(pEl,
2984  ManglerRotor,
2985  ManglerRotor(uLabel, pDO, pCraft, rrot,
2986  pRotor, pGround,
2987  ppres, dOR, dR, pdW,
2988  iMaxIter, dTolerance, dEta,
2989  dCH, dCFF,
2990  fOut));
2991  break;
2992 
2993  case DYNAMICINFLOW:
2994  DEBUGCOUT("Dynamic inflow" << std::endl);
2995 
2996  SAFENEWWITHCONSTRUCTOR(pEl,
2997  DynamicInflowRotor,
2998  DynamicInflowRotor(uLabel, pDO,
2999  pCraft, rrot, pRotor,
3000  pGround, ppres,
3001  dOR, dR,
3002  iMaxIter, dTolerance, dEta,
3003  dCH, dCFF,
3004  dVConst, dVSine, dVCosine,
3005  fOut));
3006  break;
3007 
3008  default:
3009  ASSERTMSG(0, "You shouldn't have reached this point");
3010  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3011  }
3012  break;
3013  }
3014 
3015  case PETERS_HE: {
3016  doublereal dOR = HP.GetReal();
3017  DEBUGCOUT("Reference rotation speed: " << dOR << std::endl);
3018  if (dOR <= 0.) {
3019  silent_cerr("Rotor(" << uLabel << "): "
3020  "invalid reference speed " << dOR
3021  << " at line " << HP.GetLineData()
3022  << std::endl);
3023  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3024  }
3025 
3026  doublereal dR = HP.GetReal();
3027  DEBUGCOUT("Radius: " << dR << std::endl);
3028  if (dR <= 0.) {
3029  silent_cerr("Rotor(" << uLabel << "): "
3030  "invalid radius " << dR
3031  << " at line " << HP.GetLineData()
3032  << std::endl);
3033  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3034  }
3035 
3036  // optional parameters
3037  bool bGotInitialValues(false);
3038  doublereal dVConst = 0.;
3039  doublereal dVSine = 0.;
3040  doublereal dVCosine = 0.;
3041  const StructNode *pGround = 0;
3042  unsigned iMaxIter = unsigned(-1);
3043  doublereal dTolerance = std::numeric_limits<double>::max();
3044  doublereal dEta = -1.;
3045  doublereal dCH = -1.;
3046  doublereal dCFF = -1.;
3047 
3048  while (HP.IsArg()) {
3049  if (HP.IsKeyWord("ground")) {
3050  /*
3051  * ground node for ground effect modeling
3052  */
3053  if (pGround != 0) {
3054  silent_cerr("Rotor(" << uLabel << "): "
3055  "providing another \"ground\" node "
3056  "at line " << HP.GetLineData()
3057  << std::endl);
3058  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3059  }
3060 
3061  /* ground node */
3062  pGround = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
3063 
3064  } else if (HP.IsKeyWord("initial" "value")) {
3065  if (bGotInitialValues) {
3066  silent_cerr("Rotor(" << uLabel << "): "
3067  "providing \"initial value\" another time "
3068  "at line " << HP.GetLineData()
3069  << std::endl);
3070  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3071  }
3072 
3073  dVConst = HP.GetReal();
3074  dVSine = HP.GetReal();
3075  dVCosine = HP.GetReal();
3076 
3077  bGotInitialValues = true;
3078 
3079  } else if (HP.IsKeyWord("max" "iterations")) {
3080  if (iMaxIter != unsigned(-1)) {
3081  silent_cerr("Rotor(" << uLabel << "): "
3082  "providing another \"max iterations\" value "
3083  "at line " << HP.GetLineData()
3084  << std::endl);
3085  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3086  }
3087 
3088  /* max iterations when computing reference inflow velocity;
3089  * after iMaxIter iterations, the current value is accepted
3090  * regardless of convergence; thus, 1 reproduces original
3091  * behavior */
3092  int i = HP.GetInt();
3093  if (i <= 0) {
3094  silent_cerr("illegal max iterations "
3095  << i << " for Rotor(" << uLabel << ")");
3096  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3097  }
3098  iMaxIter = i;
3099 
3100  } else if (HP.IsKeyWord("tolerance")) {
3101  if (dTolerance != std::numeric_limits<double>::max()) {
3102  silent_cerr("Rotor(" << uLabel << "): "
3103  "providing another \"tolerance\" value "
3104  "at line " << HP.GetLineData()
3105  << std::endl);
3106  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3107  }
3108 
3109  /* tolerance when computing reference inflow velocity;
3110  * when the difference in inflow velocity between two
3111  * iterations is less than tolerance in module, the
3112  * cycle breaks */
3113 
3114  dTolerance = HP.GetReal();
3115  if (dTolerance <= 0.) {
3116  silent_cerr("illegal tolerance "
3117  << dTolerance << " for Rotor(" << uLabel << ")");
3118  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3119  }
3120 
3121  } else if (HP.IsKeyWord("eta")) {
3122  if (dEta != -1.) {
3123  silent_cerr("Rotor(" << uLabel << "): "
3124  "providing another \"eta\" relaxation factor "
3125  "at line " << HP.GetLineData()
3126  << std::endl);
3127  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3128  }
3129 
3130  /* increment factor when computing reference inflow velocity;
3131  * only a fraction dEta of the difference between two iterations
3132  * is applied */
3133 
3134  dEta = HP.GetReal();
3135  if (dEta <= 0.) {
3136  silent_cerr("illegal eta "
3137  << dEta << " for Rotor(" << uLabel << ")");
3138  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3139  }
3140 
3141  } else if (HP.IsKeyWord("correction")) {
3142  if (dCH != -1.) {
3143  silent_cerr("Rotor(" << uLabel << "): "
3144  "providing another \"correction\" factor "
3145  "at line " << HP.GetLineData()
3146  << std::endl);
3147  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3148  }
3149 
3150  /* Legge la correzione della velocita' indotta */
3151  dCH = HP.GetReal();
3152  DEBUGCOUT("Hover correction: " << dCH << std::endl);
3153  if (dCH <= 0.) {
3154  silent_cerr("Rotor(" << uLabel << "): "
3155  "illegal null or negative hover inflow correction "
3156  "at line " << HP.GetLineData()
3157  << std::endl);
3158  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3159  }
3160 
3161  dCFF = HP.GetReal();
3162  DEBUGCOUT("Forward-flight correction: " << dCFF << std::endl);
3163  if (dCFF <= 0.) {
3164  silent_cerr("Rotor(" << uLabel << "): "
3165  "illegal null or negative forward-flight inflow correction "
3166  "at line " << HP.GetLineData()
3167  << std::endl);
3168  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3169  }
3170 
3171  } else {
3172  break;
3173  }
3174  }
3175 
3176  ppres = ReadResSets(pDM, HP);
3177 
3178  flag fOut = pDM->fReadOutput(HP, Elem::INDUCEDVELOCITY);
3179 
3180  if (iMaxIter == unsigned(-1)) {
3181  iMaxIter = 1;
3182 
3183  } else {
3184  if (dTolerance == std::numeric_limits<double>::max()) {
3185  silent_cerr("Rotor(" << uLabel << "): "
3186  "warning, \"max iterations\" is meaningless with default tolerance"
3187  << std::endl);
3188  }
3189  }
3190 
3191  if (dEta == -1.) {
3192  dEta = 1.;
3193  }
3194 
3195  if (dCH == -1.) {
3196  dCH = 1.;
3197  dCFF = 1.;
3198  }
3199 
3200  // create element
3201  DEBUGCOUT("Dynamic inflow" << std::endl);
3202 
3203  SAFENEWWITHCONSTRUCTOR(pEl,
3204  PetersHeRotor,
3205  PetersHeRotor(uLabel, pDO,
3206  pCraft, rrot, pRotor,
3207  pGround, ppres,
3208  dOR, dR,
3209  iMaxIter, dTolerance, dEta,
3210  dCH, dCFF,
3211  dVConst, dVSine, dVCosine,
3212  fOut));
3213  break;
3214  }
3215 
3216  default:
3217  silent_cerr("Rotor(" << uLabel << "): "
3218  "unknown induced velocity type at line "
3219  << HP.GetLineData() << std::endl);
3220  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3221  }
3222 
3223  /* Se non c'e' il punto e virgola finale */
3224  if (HP.IsArg()) {
3225  silent_cerr("Rotor(" << uLabel << "): "
3226  "semicolon expected at line "
3227  << HP.GetLineData() << std::endl);
3228  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3229  }
3230 
3231  ASSERT(pEl != 0);
3232  return pEl;
3233 } /* End of DataManager::ReadRotor() */
3234 
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
GlauertRotor::gtype
Type gtype
Definition: rotor.h:395
InducedVelocity::Res
ExternResForces Res
Definition: indvel.h:114
AirPropOwner::dGetAirDensity
virtual doublereal dGetAirDensity(const Vec3 &X) const
Definition: aerodyn.cc:736
PetersHeRotor::SetInitialValue
virtual void SetInitialValue(VectorHandler &X)
Definition: rotor.cc:2485
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
DynamicInflowRotor::AssJac
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1755
ManglerRotor::GetInducedVelocity
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:1460
Rotor::dHoverCorrection
doublereal dHoverCorrection
Definition: rotor.h:70
Rotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:535
InducedVelocity::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: indvel.cc:170
M_PI
#define M_PI
Definition: gradienttest.cc:67
Zero3
const Vec3 Zero3(0., 0., 0.)
Rotor::dArea
doublereal dArea
Definition: rotor.h:53
Rotor::dForwardFlightCorrection
doublereal dForwardFlightCorrection
Definition: rotor.h:72
Rotor::dVelocity
doublereal dVelocity
Definition: rotor.h:86
PetersHeRotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:2511
flag
long int flag
Definition: mbdyn.h:43
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
UniformRotor2::~UniformRotor2
virtual ~UniformRotor2(void)
Definition: rotor.cc:923
PetersHeRotor::dL33
doublereal dL33
Definition: rotor.h:678
InducedVelocity::ResetForce
virtual void ResetForce(void)
Definition: indvel.cc:276
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
UniformRotor::~UniformRotor
virtual ~UniformRotor(void)
Definition: rotor.cc:774
Rotor::dGetPsi
virtual doublereal dGetPsi(const Vec3 &X) const
Definition: rotor.cc:286
Rotor::dCosAlphad
doublereal dCosAlphad
Definition: rotor.h:81
UniformRotor::UniformRotor
UniformRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:694
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
GlauertRotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:1112
dM22
static const doublereal dM22
Definition: rotor.cc:1532
grad::sin
GradientExpression< UnaryExpr< FuncSin, Expr > > sin(const GradientExpression< Expr > &u)
Definition: gradient.h:2977
InducedVelocity::ppRes
ResForceSet ** ppRes
Definition: indvel.h:116
Rotor::iNumSteps
int iNumSteps
Definition: rotor.h:92
HighParser::GetInt
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
Rotor::RRot3
Vec3 RRot3
Definition: rotor.h:77
UniformRotor2::AddSectionalForce
virtual void AddSectionalForce(Elem::Type type, const Elem *pEl, unsigned uPnt, const Vec3 &F, const Vec3 &M, doublereal dW, const Vec3 &X, const Mat3x3 &R, const Vec3 &V, const Vec3 &W)
Definition: rotor.cc:936
NoRotor::NoRotor
NoRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:547
UniformRotor::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:852
ResForces::Force
virtual const Vec3 & Force(void) const
Definition: resforces.cc:103
PetersHeRotor::dVSine
doublereal dVSine
Definition: rotor.h:671
PetersHeRotor::dL31
doublereal dL31
Definition: rotor.h:677
UniformRotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:842
ManglerRotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:1417
Rotor::iMaxIter
unsigned int iMaxIter
Definition: rotor.h:59
SAFEDELETEARR
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
Rotor::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: rotor.cc:152
DynamicInflowRotor::dVCosine
doublereal dVCosine
Definition: rotor.h:547
GlauertRotor::GlauertRotor
GlauertRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:976
ASSERTMSG
#define ASSERTMSG(expr, msg)
Definition: myassert.h:219
ExternResForces::Pole
const Vec3 & Pole(void) const
Definition: resforces.cc:145
Eye3
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
DynamicInflowRotor::dL22
doublereal dL22
Definition: rotor.h:551
Rotor::RRot
Mat3x3 RRot
Definition: rotor.h:76
DynamicInflowRotor::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1785
GlauertRotor::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:1123
MYSLEEP
#define MYSLEEP(t)
Definition: mysleep.h:49
SparseSubMatrixHandler::ResizeReset
void ResizeReset(integer iNewRow, integer iNewCol)
Definition: submat.cc:1084
NoRotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:645
PetersHeRotor::dVConst
doublereal dVConst
Definition: rotor.h:670
Rotor::Init
virtual void Init(const StructNode *pC, const Mat3x3 &rrot, const StructNode *pR, const StructNode *pG, ResForceSet **ppres, const doublereal &dR, unsigned int iMaxIt, const doublereal &dTol, const doublereal &dE, flag fOut)
Definition: rotor.cc:112
PetersHeRotor::dL22
doublereal dL22
Definition: rotor.h:676
NO_OP
#define NO_OP
Definition: myassert.h:74
Elem::GetElemType
virtual Elem::Type GetElemType(void) const =0
ManglerRotor::ManglerRotor
ManglerRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:1261
Rotor::dSinAlphad
doublereal dSinAlphad
Definition: rotor.h:80
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
ResForces::Moment
virtual const Vec3 & Moment(void) const
Definition: resforces.cc:109
PetersHeRotor::Init
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, const StructNode *pGround, ResForceSet **ppres, const doublereal &dOR, const doublereal &dR, unsigned int iMaxIt, const doublereal &dTol, const doublereal &dE, const doublereal &dCH, const doublereal &dCFF, const doublereal &dVConstTmp, const doublereal &dVSineTmp, const doublereal &dVCosineTmp, flag fOut)
Definition: rotor.cc:2093
Rotor::Weight
DriveOwner Weight
Definition: rotor.h:65
Rotor::~Rotor
virtual ~Rotor(void)
Definition: rotor.cc:146
PetersHeRotor::AssJac
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:2273
DynamicInflowRotor::dVConst
doublereal dVConst
Definition: rotor.h:545
Mat3x3::GetVec
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
UniformRotor::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:785
DynamicInflowRotor::GetInducedVelocity
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:2037
NoRotor::~NoRotor
virtual ~NoRotor(void)
Definition: rotor.cc:596
DriveCaller::Restart
virtual std::ostream & Restart(std::ostream &out) const =0
DynamicInflowRotor::SetInitialValue
virtual void SetInitialValue(VectorHandler &X)
Definition: rotor.cc:1967
ManglerRotor::Init
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, const StructNode *pGround, ResForceSet **ppres, const doublereal &dOR, const doublereal &dR, DriveCaller *pdW, unsigned int iMaxIt, const doublereal &dTol, const doublereal &dE, const doublereal &dCH, const doublereal &dCFF, flag fOut)
Definition: rotor.cc:1291
DynamicInflowRotor::DynamicInflowRotor
DynamicInflowRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:1537
NoRotor
Definition: rotor.h:209
GlauertRotor::~GlauertRotor
virtual ~GlauertRotor(void)
Definition: rotor.cc:1067
SubVectorHandler::PutRowIndex
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
Zero3x3
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
InducedVelocity::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: indvel.cc:252
NoRotor::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:652
UniformRotor2::bSectionalForces
virtual bool bSectionalForces(void) const
Definition: rotor.cc:929
dM33
static const doublereal dM33
Definition: rotor.cc:1533
Rotor::dUMeanPrev
doublereal dUMeanPrev
Definition: rotor.h:56
PetersHeRotor::~PetersHeRotor
virtual ~PetersHeRotor(void)
Definition: rotor.cc:2158
SparseSubMatrixHandler::PutItem
void PutItem(integer iSubIt, integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:997
Rotor::bUMeanRefConverged
bool bUMeanRefConverged
Definition: rotor.h:63
ResForces::AddForces
void AddForces(const Vec3 &f, const Vec3 &c, const Vec3 &x)
Definition: resforces.cc:77
PetersHeRotor::PetersHeRotor
PetersHeRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:2055
Rotor::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: rotor.h:146
dVTipTreshold
static const doublereal dVTipTreshold
Definition: rotor.cc:55
Rotor::Rotor
Rotor(unsigned int uL, const DofOwner *pDO)
Definition: rotor.cc:59
Rotor::dOmega
doublereal dOmega
Definition: rotor.h:87
PetersHeRotor::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:2303
SAFENEW
#define SAFENEW(pnt, item)
Definition: mynewmem.h:695
PetersHeRotor::dL11
doublereal dL11
Definition: rotor.h:674
GlauertRotor::GetInducedVelocity
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:1158
DynamicInflowRotor::Output
virtual void Output(OutputHandler &OH) const
Definition: rotor.cc:1651
VariableSubMatrixHandler::SetNullMatrix
void SetNullMatrix(void)
Definition: submat.h:1159
false
Definition: mbdyn.h:76
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
DynamicInflowRotor::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:1993
PetersHeRotor::GetInducedVelocity
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:2555
Rotor::dPsi0
doublereal dPsi0
Definition: rotor.h:79
copysign
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
UniformRotor::Init
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, const StructNode *pGround, ResForceSet **ppres, const doublereal &dOR, const doublereal &dR, DriveCaller *pdW, unsigned int iMaxIt, const doublereal &dTol, const doublereal &dE, const doublereal &dCH, const doublereal &dCFF, flag fOut)
Definition: rotor.cc:725
Rotor::VCraft
Vec3 VCraft
Definition: rotor.h:78
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
Rotor::dUMeanRef
doublereal dUMeanRef
Definition: rotor.h:55
NoRotor::GetInducedVelocity
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:683
Rotor::dMu
doublereal dMu
Definition: rotor.h:82
StructNode::GetWCurr
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
Rotor::GetPos
virtual void GetPos(const Vec3 &X, doublereal &dr, doublereal &dp) const
Definition: rotor.cc:304
ExternResForces::PutPole
void PutPole(const Vec3 &x)
Definition: resforces.cc:139
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
DynamicInflowRotor::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:2003
ReadRotor
Elem * ReadRotor(DataManager *pDM, MBDynParser &HP, const DofOwner *pDO, unsigned int uLabel)
Definition: rotor.cc:2573
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
KeyWords
KeyWords
Definition: dataman4.cc:94
VectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
GlauertRotor::Init
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, const StructNode *pGround, ResForceSet **ppres, const doublereal &dOR, const doublereal &dR, DriveCaller *pdW, unsigned int iMaxIt, const doublereal &dTol, const doublereal &dE, const doublereal &dCH, const doublereal &dCFF, flag fOut)
Definition: rotor.cc:1010
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
DriveOwner::pGetDriveCaller
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
ReadResSets
ResForceSet ** ReadResSets(DataManager *pDM, MBDynParser &HP)
Definition: resforces.cc:263
PetersHeRotor::dVCosine
doublereal dVCosine
Definition: rotor.h:672
PetersHeRotor::Output
virtual void Output(OutputHandler &OH) const
Definition: rotor.cc:2169
DynamicInflowRotor::dL11
doublereal dL11
Definition: rotor.h:549
DynamicInflowRotor::dVSine
doublereal dVSine
Definition: rotor.h:546
Rotor::dVTipRef
doublereal dVTipRef
Definition: rotor.h:52
Rotor::RRotTranspose
Mat3x3 RRotTranspose
Definition: rotor.h:75
Rotor::dGetPos
virtual doublereal dGetPos(const Vec3 &X) const
Definition: rotor.cc:296
grad::Norm
VectorExpression< VectorExpr, N_rows >::ScalarType Norm(const VectorExpression< VectorExpr, N_rows > &u)
Definition: matvec.h:3163
Mat3x3::Transpose
Mat3x3 Transpose(void) const
Definition: matvec3.h:816
Mat3x3::pGetMat
const doublereal * pGetMat(void) const
Definition: matvec3.h:743
ResForces::AddForce
void AddForce(const Vec3 &f)
Definition: resforces.cc:58
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
Rotor::dEta
doublereal dEta
Definition: rotor.h:62
Elem
Definition: elem.h:75
DynamicInflowRotor::dL13
doublereal dL13
Definition: rotor.h:550
Elem::Type
Type
Definition: elem.h:91
Rotor::dOmegaRef
doublereal dOmegaRef
Definition: rotor.h:49
HighParser::GetWord
virtual int GetWord(void)
Definition: parser.cc:1083
OutputHandler::Rotors
std::ostream & Rotors(void) const
Definition: output.h:464
Vec3::pGetVec
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
DynamicInflowRotor::~DynamicInflowRotor
virtual ~DynamicInflowRotor(void)
Definition: rotor.cc:1640
ToBeOutput::fToBeOutput
virtual flag fToBeOutput(void) const
Definition: output.cc:884
DriveOwner::Set
void Set(const DriveCaller *pDC)
Definition: drive.cc:647
DriveOwner::dGet
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
Rotor::InitParam
virtual void InitParam(bool bComputeMeanInducedVelocity=true)
Definition: rotor.cc:323
Rotor::iCurrIter
unsigned int iCurrIter
Definition: rotor.h:60
Rotor::dWeight
doublereal dWeight
Definition: rotor.h:68
Rotor::dLambda
doublereal dLambda
Definition: rotor.h:83
PetersHeRotor::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:2521
NoRotor::Init
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, ResForceSet **ppres, const doublereal &dR, flag fOut)
Definition: rotor.cc:570
StructDispNode::GetVCurr
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
UniformRotor2::UniformRotor2
UniformRotor2(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:894
Rotor::dRadius
doublereal dRadius
Definition: rotor.h:51
SAFENEWARR
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
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
ManglerRotor::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1358
MBDynParser::GetDriveCaller
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
Rotor::pRotor
const StructNode * pRotor
Definition: rotor.h:46
DynamicInflowRotor::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: rotor.cc:1975
PetersHeRotor::dL13
doublereal dL13
Definition: rotor.h:675
DynamicInflowRotor::dL33
doublereal dL33
Definition: rotor.h:553
GlauertRotor::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1079
VariableSubMatrixHandler::SetSparse
SparseSubMatrixHandler & SetSparse(void)
Definition: submat.h:1178
ToBeOutput::SetOutputFlag
virtual void SetOutputFlag(flag f=flag(1))
Definition: output.cc:896
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
grad::atan2
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
DynamicInflowRotor::dL31
doublereal dL31
Definition: rotor.h:552
doublereal
double doublereal
Definition: colamd.c:52
Rotor::Output
virtual void Output(OutputHandler &OH) const
Definition: rotor.cc:183
integer
long int integer
Definition: colamd.c:51
dM11
static const doublereal dM11
Definition: rotor.cc:1529
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
Rotor::pGround
const StructNode * pGround
Definition: rotor.h:47
PetersHeRotor::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: rotor.cc:2493
ManglerRotor::AddForce
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:1427
UniformRotor::GetInducedVelocity
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:884
grad::tan
GradientExpression< UnaryExpr< FuncTan, Expr > > tan(const GradientExpression< Expr > &u)
Definition: gradient.h:2979
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Rotor::dUMean
doublereal dUMean
Definition: rotor.h:54
DataManager::ReadNode
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
Rotor::dTolerance
doublereal dTolerance
Definition: rotor.h:61
Rotor::dChi
doublereal dChi
Definition: rotor.h:84
DynamicInflowRotor::Init
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, const StructNode *pGround, ResForceSet **ppres, const doublereal &dOR, const doublereal &dR, unsigned int iMaxIt, const doublereal &dTol, const doublereal &dE, const doublereal &dCH, const doublereal &dCFF, const doublereal &dVConstTmp, const doublereal &dVSineTmp, const doublereal &dVCosineTmp, flag fOut)
Definition: rotor.cc:1575
VectorHandler::Resize
virtual void Resize(integer iNewSize)=0
R
Mat3x3 R
Definition: mbdyn_util_oct.cc:74
InducedVelocity::pCraft
const StructNode * pCraft
Definition: indvel.h:111
Rotor
Definition: rotor.h:43
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
NoRotor::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:607
ManglerRotor::~ManglerRotor
virtual ~ManglerRotor(void)
Definition: rotor.cc:1346
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056