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 {
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 
147 {
148  NO_OP;
149 }
150 
151 void
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 
180 }
181 
182 void
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  */
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 */
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 */
332 
333  /* Posizione del rotore */
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);
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
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) {
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);
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) {
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 
597 {
598 #ifdef USE_MPI
599  SAFEDELETEARR(pBlockLenght);
600  SAFEDELETEARR(pDispl);
601  SAFEDELETE(pIndVelDataType);
602 #endif /* USE_MPI */
603 }
604 
605 /* assemblaggio residuo */
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
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
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;
748  dArea = M_PI*dRadius*dRadius;
749  Weight.Set(pdW);
750  dHoverCorrection = dCH;
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 
775 {
776 #ifdef USE_MPI
777  SAFEDELETEARR(pBlockLenght);
778  SAFEDELETEARR(pDispl);
779  SAFEDELETE(pIndVelDataType);
780 #endif /* USE_MPI */
781 }
782 
783 /* assemblaggio residuo */
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 */
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 << ", ",
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
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 
924 {
925  NO_OP;
926 }
927 
928 bool
930 {
931  return true;
932 }
933 
934 /* Somma alla trazione il contributo di forza di un elemento generico */
935 void
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
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;
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 
1068 {
1069 #ifdef USE_MPI
1070  SAFEDELETEARR(pBlockLenght);
1071  SAFEDELETEARR(pDispl);
1072  SAFEDELETE(pIndVelDataType);
1073 #endif /* USE_MPI */
1074 }
1075 
1076 
1077 /* assemblaggio residuo */
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 << ", ",
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
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:
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
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;
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 
1347 {
1348 #ifdef USE_MPI
1349  SAFEDELETEARR(pBlockLenght);
1350  SAFEDELETEARR(pDispl);
1351  SAFEDELETE(pIndVelDataType);
1352 #endif /* USE_MPI */
1353 }
1354 
1355 
1356 /* assemblaggio residuo */
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 << ", ",
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
1461  unsigned uLabel, unsigned uPnt, const Vec3& X) const
1462 {
1463  if (dUMeanPrev == 0.) {
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 
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
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;
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 
1641 {
1642 #ifdef USE_MPI
1643  SAFEDELETEARR(pBlockLenght);
1644  SAFEDELETEARR(pDispl);
1645  SAFEDELETE(pIndVelDataType);
1646 #endif /* USE_MPI */
1647 }
1648 
1649 
1650 void
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 */
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 */
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 
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
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  */
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
1968 {
1969  NO_OP;
1970 }
1971 
1972 
1973 /* Relativo ai ...WithDofs */
1974 void
1976  VectorHandler& X, VectorHandler& XP,
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
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
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
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;
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 
2159 {
2160 #ifdef USE_MPI
2161  SAFEDELETEARR(pBlockLenght);
2162  SAFEDELETEARR(pDispl);
2163  SAFEDELETE(pIndVelDataType);
2164 #endif /* USE_MPI */
2165 }
2166 
2167 
2168 void
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 */
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 */
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 
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
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  */
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
2486 {
2487  NO_OP;
2488 }
2489 
2490 
2491 /* Relativo ai ...WithDofs */
2492 void
2494  VectorHandler& X, VectorHandler& XP,
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
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 *
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);
2655  }
2656  }
2657 
2658  ppres = ReadResSets(pDM, HP);
2659 
2660  flag fOut = pDM->fReadOutput(HP, Elem::INDUCEDVELOCITY);
2661 
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: {
2674  if (InducedType == GLAUERT && HP.IsKeyWord("type")) {
2675  if (HP.IsKeyWord("glauert")) {
2676  gtype = GlauertRotor::GLAUERT;
2677 
2678  } else if (HP.IsKeyWord("coleman")) {
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")) {
2689 
2690  } else if (HP.IsKeyWord("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);
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);
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);
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);
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);
2764  }
2765 
2766  if (bGotInitialValues) {
2767  silent_cerr("Rotor(" << uLabel << "): "
2768  "providing \"initial value\" another time "
2769  "at line " << HP.GetLineData()
2770  << std::endl);
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);
2787  }
2788 
2789  if (pdW != 0) {
2790  silent_cerr("Rotor(" << uLabel << "): "
2791  "providing another \"delay\" driver "
2792  "at line " << HP.GetLineData()
2793  << std::endl);
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);
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 << ")");
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);
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 << ")");
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);
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 << ")");
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);
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);
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);
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);
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);
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);
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 
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 
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");
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);
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);
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);
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);
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);
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 << ")");
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);
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 << ")");
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);
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 << ")");
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);
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);
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);
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 
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);
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);
3229  }
3230 
3231  ASSERT(pEl != 0);
3232  return pEl;
3233 } /* End of DataManager::ReadRotor() */
3234 
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
Type gtype
Definition: rotor.h:395
ExternResForces Res
Definition: indvel.h:114
virtual doublereal dGetAirDensity(const Vec3 &X) const
Definition: aerodyn.cc:736
virtual void SetInitialValue(VectorHandler &X)
Definition: rotor.cc:2485
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1755
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:1460
doublereal dHoverCorrection
Definition: rotor.h:70
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:535
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: indvel.cc:170
#define M_PI
Definition: gradienttest.cc:67
const Vec3 Zero3(0., 0., 0.)
doublereal dArea
Definition: rotor.h:53
doublereal dForwardFlightCorrection
Definition: rotor.h:72
doublereal dVelocity
Definition: rotor.h:86
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:2511
long int flag
Definition: mbdyn.h:43
virtual bool bToBeOutput(void) const
Definition: output.cc:890
virtual ~UniformRotor2(void)
Definition: rotor.cc:923
doublereal dL33
Definition: rotor.h:678
virtual void ResetForce(void)
Definition: indvel.cc:276
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
virtual ~UniformRotor(void)
Definition: rotor.cc:774
virtual doublereal dGetPsi(const Vec3 &X) const
Definition: rotor.cc:286
doublereal dCosAlphad
Definition: rotor.h:81
UniformRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:694
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:1112
static const doublereal dM22
Definition: rotor.cc:1532
GradientExpression< UnaryExpr< FuncSin, Expr > > sin(const GradientExpression< Expr > &u)
Definition: gradient.h:2977
ResForceSet ** ppRes
Definition: indvel.h:116
int iNumSteps
Definition: rotor.h:92
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
Vec3 RRot3
Definition: rotor.h:77
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(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:547
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:852
virtual const Vec3 & Force(void) const
Definition: resforces.cc:103
doublereal dVSine
Definition: rotor.h:671
doublereal dL31
Definition: rotor.h:677
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:842
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:1417
unsigned int iMaxIter
Definition: rotor.h:59
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: rotor.cc:152
doublereal dVCosine
Definition: rotor.h:547
GlauertRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:976
#define ASSERTMSG(expr, msg)
Definition: myassert.h:219
const Vec3 & Pole(void) const
Definition: resforces.cc:145
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
doublereal dL22
Definition: rotor.h:551
Mat3x3 RRot
Definition: rotor.h:76
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1785
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:1123
#define MYSLEEP(t)
Definition: mysleep.h:49
void ResizeReset(integer iNewRow, integer iNewCol)
Definition: submat.cc:1084
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:645
doublereal dVConst
Definition: rotor.h:670
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
doublereal dL22
Definition: rotor.h:676
#define NO_OP
Definition: myassert.h:74
virtual Elem::Type GetElemType(void) const =0
ManglerRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:1261
doublereal dSinAlphad
Definition: rotor.h:80
std::vector< Hint * > Hints
Definition: simentity.h:89
virtual const Vec3 & Moment(void) const
Definition: resforces.cc:109
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
DriveOwner Weight
Definition: rotor.h:65
virtual ~Rotor(void)
Definition: rotor.cc:146
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:2273
doublereal dVConst
Definition: rotor.h:545
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:785
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:2037
virtual ~NoRotor(void)
Definition: rotor.cc:596
virtual std::ostream & Restart(std::ostream &out) const =0
virtual void SetInitialValue(VectorHandler &X)
Definition: rotor.cc:1967
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(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:1537
Definition: rotor.h:209
virtual ~GlauertRotor(void)
Definition: rotor.cc:1067
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: indvel.cc:252
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:652
virtual bool bSectionalForces(void) const
Definition: rotor.cc:929
static const doublereal dM33
Definition: rotor.cc:1533
doublereal dUMeanPrev
Definition: rotor.h:56
virtual ~PetersHeRotor(void)
Definition: rotor.cc:2158
void PutItem(integer iSubIt, integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:997
bool bUMeanRefConverged
Definition: rotor.h:63
void AddForces(const Vec3 &f, const Vec3 &c, const Vec3 &x)
Definition: resforces.cc:77
PetersHeRotor(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:2055
virtual const Vec3 & GetXCurr(void) const
Definition: rotor.h:146
static const doublereal dVTipTreshold
Definition: rotor.cc:55
Rotor(unsigned int uL, const DofOwner *pDO)
Definition: rotor.cc:59
doublereal dOmega
Definition: rotor.h:87
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:2303
#define SAFENEW(pnt, item)
Definition: mynewmem.h:695
doublereal dL11
Definition: rotor.h:674
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:1158
virtual void Output(OutputHandler &OH) const
Definition: rotor.cc:1651
void SetNullMatrix(void)
Definition: submat.h:1159
Definition: mbdyn.h:76
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
virtual std::ostream & Restart(std::ostream &out) const
Definition: rotor.cc:1993
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:2555
doublereal dPsi0
Definition: rotor.h:79
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
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
Vec3 VCraft
Definition: rotor.h:78
#define DEBUGCOUT(msg)
Definition: myassert.h:232
doublereal dUMeanRef
Definition: rotor.h:55
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:683
doublereal dMu
Definition: rotor.h:82
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
virtual void GetPos(const Vec3 &X, doublereal &dr, doublereal &dp) const
Definition: rotor.cc:304
void PutPole(const Vec3 &x)
Definition: resforces.cc:139
#define ASSERT(expression)
Definition: colamd.c:977
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:2003
Elem * ReadRotor(DataManager *pDM, MBDynParser &HP, const DofOwner *pDO, unsigned int uLabel)
Definition: rotor.cc:2573
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
KeyWords
Definition: dataman4.cc:94
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
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
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
ResForceSet ** ReadResSets(DataManager *pDM, MBDynParser &HP)
Definition: resforces.cc:263
doublereal dVCosine
Definition: rotor.h:672
virtual void Output(OutputHandler &OH) const
Definition: rotor.cc:2169
doublereal dL11
Definition: rotor.h:549
doublereal dVSine
Definition: rotor.h:546
doublereal dVTipRef
Definition: rotor.h:52
Mat3x3 RRotTranspose
Definition: rotor.h:75
virtual doublereal dGetPos(const Vec3 &X) const
Definition: rotor.cc:296
VectorExpression< VectorExpr, N_rows >::ScalarType Norm(const VectorExpression< VectorExpr, N_rows > &u)
Definition: matvec.h:3163
Mat3x3 Transpose(void) const
Definition: matvec3.h:816
const doublereal * pGetMat(void) const
Definition: matvec3.h:743
void AddForce(const Vec3 &f)
Definition: resforces.cc:58
virtual bool IsArg(void)
Definition: parser.cc:807
doublereal dEta
Definition: rotor.h:62
Definition: elem.h:75
doublereal dL13
Definition: rotor.h:550
Type
Definition: elem.h:91
doublereal dOmegaRef
Definition: rotor.h:49
virtual int GetWord(void)
Definition: parser.cc:1083
std::ostream & Rotors(void) const
Definition: output.h:464
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
virtual ~DynamicInflowRotor(void)
Definition: rotor.cc:1640
virtual flag fToBeOutput(void) const
Definition: output.cc:884
void Set(const DriveCaller *pDC)
Definition: drive.cc:647
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
virtual void InitParam(bool bComputeMeanInducedVelocity=true)
Definition: rotor.cc:323
unsigned int iCurrIter
Definition: rotor.h:60
doublereal dWeight
Definition: rotor.h:68
doublereal dLambda
Definition: rotor.h:83
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:2521
virtual void Init(const StructNode *pCraft, const Mat3x3 &rrot, const StructNode *pRotor, ResForceSet **ppres, const doublereal &dR, flag fOut)
Definition: rotor.cc:570
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
UniformRotor2(unsigned int uLabel, const DofOwner *pDO)
Definition: rotor.cc:894
doublereal dRadius
Definition: rotor.h:51
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
virtual flag fGetAirVelocity(Vec3 &Velocity, const Vec3 &X) const
Definition: aerodyn.cc:725
GradientExpression< UnaryExpr< FuncCos, Expr > > cos(const GradientExpression< Expr > &u)
Definition: gradient.h:2978
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1358
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
const StructNode * pRotor
Definition: rotor.h:46
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: rotor.cc:1975
doublereal dL13
Definition: rotor.h:675
doublereal dL33
Definition: rotor.h:553
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:1079
SparseSubMatrixHandler & SetSparse(void)
Definition: submat.h:1178
virtual void SetOutputFlag(flag f=flag(1))
Definition: output.cc:896
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
doublereal dL31
Definition: rotor.h:552
double doublereal
Definition: colamd.c:52
virtual void Output(OutputHandler &OH) const
Definition: rotor.cc:183
long int integer
Definition: colamd.c:51
static const doublereal dM11
Definition: rotor.cc:1529
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
const StructNode * pGround
Definition: rotor.h:47
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: rotor.cc:2493
virtual void AddForce(const Elem *pEl, const StructNode *pNode, const Vec3 &F, const Vec3 &M, const Vec3 &X)
Definition: rotor.cc:1427
virtual Vec3 GetInducedVelocity(Elem::Type type, unsigned uLabel, unsigned uPnt, const Vec3 &X) const
Definition: rotor.cc:884
GradientExpression< UnaryExpr< FuncTan, Expr > > tan(const GradientExpression< Expr > &u)
Definition: gradient.h:2979
unsigned int GetLabel(void) const
Definition: withlab.cc:62
doublereal dUMean
Definition: rotor.h:54
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
doublereal dTolerance
Definition: rotor.h:61
doublereal dChi
Definition: rotor.h:84
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
virtual void Resize(integer iNewSize)=0
Mat3x3 R
const StructNode * pCraft
Definition: indvel.h:111
Definition: rotor.h:43
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: rotor.cc:607
virtual ~ManglerRotor(void)
Definition: rotor.cc:1346
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056