MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
aeromodal.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aeromodal.cc,v 1.64 2017/01/12 14:45:58 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31  /* Elemento aerodinamico modale */
32 
33 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
34 
35 #include <cerrno>
36 
37 #include "aeromodal.h"
38 #include "dataman.h"
39 #include "Rot.hh"
40 
41 /* AerodynamicModal - begin */
42 
43 AerodynamicModal::AerodynamicModal(unsigned int uLabel,
44  const StructNode* pN,
45  const Modal* pMJ,
46  const Mat3x3& RaTmp,
47  const DofOwner* pDO,
48  doublereal Cd,
49  const int NModal,
50  const int NAero,
51  RigidF_t rgF,
52  const int Gust,
53  const doublereal Vff,
54  SpMapMatrixHandler* pAMat,
55  FullMatrixHandler* pBMat,
56  FullMatrixHandler* pCMat,
57  FullMatrixHandler* pD0Mat,
58  FullMatrixHandler* pD1Mat,
59  FullMatrixHandler* pD2Mat,
60  flag fout)
61 : Elem(uLabel, fout),
62 AerodynamicElem(uLabel, pDO, fout),
63 InitialAssemblyElem(uLabel, fout),
64 pModalNode(pN), pModalJoint(pMJ),
65 Ra(RaTmp),
66 Chord(Cd),
67 NStModes(NModal),
68 NAeroStates(NAero),
69 NGust(Gust),
70 pA(pAMat), pB(pBMat), pC(pCMat),
71 pD0(pD0Mat), pD1(pD1Mat), pD2(pD2Mat),
72 pq(0), pqPrime(0), pqSec(0),
73 pxa(0), pxaPrime(0),
74 pgs(0), pgsPrime(0),
75 gustVff(Vff), gustXi(0.707),
76 RigidF(rgF)
77 {
78  DEBUGCOUTFNAME("AerodynamicModal::AerodynamicModal");
79 
80  R0 = pModalNode->GetRCurr()*Ra;
81  P0 = R0.MulTV(pModalNode->GetXCurr());
82 
83  SAFENEWWITHCONSTRUCTOR(pq, MyVectorHandler,
84  MyVectorHandler(NStModes+RigidF));
85  SAFENEWWITHCONSTRUCTOR(pqPrime, MyVectorHandler,
86  MyVectorHandler(NStModes+RigidF));
87  SAFENEWWITHCONSTRUCTOR(pqSec, MyVectorHandler,
88  MyVectorHandler(NStModes+RigidF));
89  SAFENEWWITHCONSTRUCTOR(pxa, MyVectorHandler,
90  MyVectorHandler(NAeroStates));
91  SAFENEWWITHCONSTRUCTOR(pxaPrime, MyVectorHandler,
92  MyVectorHandler(NAeroStates));
93  SAFENEWWITHCONSTRUCTOR(pgs, MyVectorHandler,
94  MyVectorHandler(2*NGust));
95  SAFENEWWITHCONSTRUCTOR(pgsPrime, MyVectorHandler,
96  MyVectorHandler(2*NGust));
97 
98  ASSERT(pModalNode != 0);
99  ASSERT(pModalNode->GetNodeType() == Node::STRUCTURAL);
100 }
101 
102 AerodynamicModal::~AerodynamicModal(void)
103 {
104  DEBUGCOUTFNAME("AerodynamicModal::~AerodynamicModal");
105  if (pq != 0) {
106  SAFEDELETE(pq);
107  }
108  if (pqPrime != 0) {
109  SAFEDELETE(pqPrime);
110  }
111  if (pqSec != 0) {
112  SAFEDELETE(pqSec);
113  }
114  if (pxa != 0) {
115  SAFEDELETE(pxa);
116  }
117  if (pxaPrime != 0) {
118  SAFEDELETE(pxaPrime);
119  }
120  if (pA != 0) {
121  SAFEDELETE(pA);
122  }
123  if (pB != 0) {
124  SAFEDELETE(pB);
125  }
126  if (pC != 0) {
127  SAFEDELETE(pC);
128  }
129  if (pD0 != 0) {
130  SAFEDELETE(pD0);
131  }
132  if (pD1 != 0) {
133  SAFEDELETE(pD1);
134  }
135  if (pD2 != 0) {
136  SAFEDELETE(pD2);
137  }
138 }
139 
140 /* Scrive il contributo dell'elemento al file di restart */
141 std::ostream&
142 AerodynamicModal::Restart(std::ostream& out) const
143 {
144  return out << " /* aerodynamic modal: not implemented yet */" << std::endl;
145 }
146 
147 /* assemblaggio jacobiano */
148 VariableSubMatrixHandler&
149 AerodynamicModal::AssJac(VariableSubMatrixHandler& WorkMat,
150  doublereal dCoef,
151  const VectorHandler& XCurr,
152  const VectorHandler& XPrimeCurr )
153 {
154  DEBUGCOUT("Entering AerodynamicModal::AssJac()" << std::endl);
155 
156  FullSubMatrixHandler& WM = WorkMat.SetFull();
157 
158  integer iNumRows = 0;
159  integer iNumCols = 0;
160  WorkSpaceDim(&iNumRows, &iNumCols);
161  WM.ResizeReset(iNumRows, iNumCols);
162 
163  integer iModalIndex = pModalJoint->iGetModalIndex();
164  for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
165  WM.PutRowIndex(iCnt, iModalIndex+NStModes+iCnt);
166  WM.PutColIndex(iCnt, iModalIndex+iCnt);
167  WM.PutColIndex(NStModes+iCnt, iModalIndex+NStModes+iCnt);
168  }
169 
170  integer iAeroIndex = iGetFirstIndex();
171  for (unsigned int iCnt = 1; iCnt <= NAeroStates+2*NGust; iCnt++) {
172  WM.PutRowIndex(NStModes+iCnt, iAeroIndex+iCnt);
173  WM.PutColIndex(2*NStModes+iCnt, iAeroIndex+iCnt);
174  }
175 
176  /* Dati del nodo rigido */
177  const Vec3& X0(pModalNode->GetXCurr());
178  const Vec3& V0(pModalNode->GetVCurr());
179  const Mat3x3& Rn(pModalNode->GetRCurr());
180  Mat3x3 RR(Rn*Ra);
181  Vec3 Vr(V0);
182 
183  doublereal rho, vs, p, T;
184  GetAirProps(X0, rho, vs, p, T); /* p, T are not used yet */
185 
186  Vec3 VTmp(Zero3);
187  if (fGetAirVelocity(VTmp, X0)) {
188  Vr -= VTmp;
189  }
190 
191  /* velocità nel riferimento nodale aerodinamico */
192  VTmp = RR.MulTV(Vr);
193  doublereal nV = std::abs(VTmp(1));
194  /* doublereal CV=Chord/(2*nV); */
195  doublereal qd = 0.5*rho*nV*nV;
196  doublereal qd1 = 0.25*rho*nV*Chord; /* qd*CV */
197  doublereal qd2 = 0.125*rho*Chord*Chord; /* qd*CV*CV */
198 
199  /* parte deformabile :
200  *
201  * | ||aP|
202  * | || |
203  * | cKae Mae+cCae -cqC ||bP|
204  * | || |
205  * |-c(1/CV)B 0 I-c(1/CV)A ||xa|
206  *
207  * con
208  * Mae=-qd CV^2 D2
209  * Cae=-qd CV D1
210  * Kae=-qd D0
211  */
212  for (unsigned int i = 1; i <= NStModes; i++) {
213  for (unsigned int j = 1; j <= NStModes; j++) {
214  WM.IncCoef(i, j,
215  -dCoef*qd*pD0->operator()(RigidF + i, RigidF + j));
216  WM.IncCoef(i, j + NStModes,
217  -qd2*pD2->operator()(RigidF + i, RigidF + j)
218  -dCoef*qd1*pD1->operator()(RigidF + i, RigidF + j));
219  }
220  }
221 
222  for (unsigned int j = 1; j <= NAeroStates; j++) {
223  for (unsigned int i = 1; i <= NStModes; i++) {
224  WM.IncCoef(i, j + 2*NStModes,
225  -dCoef*qd*pC->operator()(RigidF + i, j));
226  WM.IncCoef(j + NStModes, i,
227  -dCoef*(2*nV/Chord)*pB->operator()(j, RigidF + i));
228  }
229  }
230 
231  for (unsigned int j = 1; j <= NAeroStates; j++) {
232  for (unsigned int i = 1; i <= NAeroStates; i++) {
233  WM.IncCoef(i + NStModes, j + 2*NStModes,
234  1.*(i == j) - dCoef*(2*nV/Chord)*pA->operator()(i, j));
235  }
236  }
237 
238  if (NGust) {
239  for (unsigned int i = 1; i <= NStModes; i++) {
240  WM.IncCoef(i, 2*NStModes + NAeroStates + 1,
241  -dCoef*qd*pD0->operator()(RigidF + i, RigidF + NStModes + 1));
242  WM.IncCoef(i, 2*NStModes + NAeroStates + 3,
243  -dCoef*qd*pD0->operator()(RigidF + i, RigidF + NStModes + 2));
244  WM.IncCoef(i, 2*NStModes + NAeroStates + 2,
245  -qd2*pD2->operator()(RigidF + i, RigidF + NStModes + 1)
246  -dCoef*qd1*pD1->operator()(RigidF + i, RigidF + NStModes + 1));
247  WM.IncCoef(i, 2*NStModes + NAeroStates + 4,
248  -qd2*pD2->operator()(RigidF + i, RigidF + NStModes + 2)
249  -dCoef*qd1*pD1->operator()(RigidF + i, RigidF + NStModes + 2));
250  }
251 
252  for (unsigned int i = 1; i <= NAeroStates; i++) {
253  WM.IncCoef(i + NStModes, 2*NStModes + NAeroStates + 1,
254  -dCoef*(2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 1));
255  WM.IncCoef(i + NStModes, 2*NStModes + NAeroStates + 3,
256  -dCoef*(2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 2));
257  }
258 
259  for (unsigned int i = 0; i < NGust; i++) {
260  WM.IncCoef(i*2 + 1 + NStModes+NAeroStates,
261  i*2 + 1 + 2*NStModes + NAeroStates, 1.);
262  WM.IncCoef(i*2 + 1 + NStModes+NAeroStates,
263  i*2 + 2 + 2*NStModes + NAeroStates, -dCoef);
264  WM.IncCoef(i*2 + 2 + NStModes+NAeroStates,
265  i*2 + 1 + 2*NStModes+NAeroStates,
266  dCoef*gustVff*gustVff);
267  WM.IncCoef(i*2 + 2 + NStModes + NAeroStates,
268  i*2 + 2 + 2*NStModes + NAeroStates,
269  1. + dCoef*2*gustXi*gustVff);
270  }
271  }
272 
273  return WorkMat;
274 }
275 
276 
277 SubVectorHandler&
278 AerodynamicModal::AssRes(SubVectorHandler& WorkVec,
279  doublereal dCoef,
280  const VectorHandler& XCurr,
281  const VectorHandler& XPrimeCurr)
282 {
283  DEBUGCOUTFNAME("AerodynamicModal::AssRes");
284  WorkVec.ResizeReset(RigidF + NStModes + NAeroStates + 2*NGust);
285 
286  if (RigidF) {
287  integer iFirstIndex = pModalNode->iGetFirstMomentumIndex();
288  for (int iCnt = 1; iCnt <= RigidF; iCnt++) {
289  WorkVec.PutRowIndex(iCnt, iFirstIndex + iCnt);
290  }
291 
292  const Vec3& X0(pModalNode->GetXCurr());
293  const Mat3x3& Rn(pModalNode->GetRCurr());
294 
295  Mat3x3 RR(Rn*Ra);
296 
297  // q
298  pq->Put(1, RR.MulTV(X0) - P0);
299  pq->Put(4, RotManip::VecRot(RR.MulMT(R0))); // nota: wrappa; se serve si può eliminare
300 
301  // dot{q}
302  const Vec3& V0(pModalNode->GetVCurr());
303  const Vec3& W0(pModalNode->GetWCurr());
304  pqPrime->Put(1, RR.MulTV(V0 + X0.Cross(W0)));
305  pqPrime->Put(4, RR.MulTV(W0));
306 
307  // ddot{q}
308  const Vec3& XPP0(pModalNode->GetXPPCurr());
309  const Vec3& WP0(pModalNode->GetWPCurr());
310  pqSec->Put(1, RR.MulTV(XPP0)); // verificare: non mancano i termini di trasporto ecc?
311  pqSec->Put(4, RR.MulTV(WP0));
312  }
313 
314  integer iModalIndex = pModalJoint->iGetModalIndex();
315 
316  for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
317  WorkVec.PutRowIndex(RigidF + iCnt, iModalIndex + NStModes + iCnt);
318  }
319 
320  integer iAeroIndex = iGetFirstIndex();
321  for (unsigned int iCnt = 1; iCnt <= NAeroStates + 2*NGust; iCnt++) {
322  WorkVec.PutRowIndex(RigidF + NStModes+iCnt, iAeroIndex + iCnt);
323  }
324 
325  /* Recupera i vettori {a} e {aP} e {aS} (deformate modali) */
326  // FIXME: get this info from modal joint?
327  for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
328  pq->PutCoef(iCnt + RigidF, XCurr(iModalIndex + iCnt));
329  pqPrime->PutCoef(iCnt + RigidF, XPrimeCurr(iModalIndex + iCnt));
330  pqSec->PutCoef(iCnt + RigidF, XPrimeCurr(iModalIndex + NStModes + iCnt));
331  }
332 
333  /* Recupera i vettori {xa} e {xaP} */
334  for (unsigned int iCnt = 1; iCnt <= NAeroStates; iCnt++) {
335  pxa->PutCoef(iCnt, XCurr(iAeroIndex + iCnt));
336  pxaPrime->PutCoef(iCnt, XPrimeCurr(iAeroIndex + iCnt));
337  }
338 
339  for (unsigned int iCnt = 1; iCnt <= 2*NGust; iCnt++) {
340  pgs->PutCoef(iCnt, XCurr(iAeroIndex + NAeroStates + iCnt));
341  pgsPrime->PutCoef(iCnt, XPrimeCurr(iAeroIndex + NAeroStates + iCnt));
342  }
343 
344  AssVec(WorkVec);
345 
346  return WorkVec;
347 }
348 
349 SubVectorHandler&
350 AerodynamicModal::InitialAssRes(SubVectorHandler& WorkVec,
351  const VectorHandler& XCurr)
352 {
353  DEBUGCOUTFNAME("AerodynamicModal::AssRes");
354  const Vec3& X0(pModalNode->GetXCurr());
355  const Mat3x3& Rn(pModalNode->GetRCurr());
356  R0 = Rn*Ra;
357  P0 = R0.MulTV(X0);
358 
359  WorkVec.ResizeReset(RigidF + NStModes + NAeroStates + 2*NGust);
360 
361  if (RigidF) {
362  integer iFirstIndex = pModalNode->iGetFirstIndex();
363  for (int iCnt = 1; iCnt <= RigidF; iCnt++) {
364  WorkVec.PutRowIndex(iCnt, iFirstIndex + RigidF + iCnt);
365  }
366 
367  const Vec3& X0(pModalNode->GetXCurr());
368  const Mat3x3& Rn(pModalNode->GetRCurr());
369  Mat3x3 RR(Rn*Ra);
370 
371  // q
372  pq->Put(1, RR.MulTV(X0) - P0);
373 
374  Vec3 g(CGR_Rot::Param, RR.MulMT(R0));
375  doublereal d(g.Norm());
376  if (d > std::numeric_limits<doublereal>::epsilon()) {
377  pq->Put(4, g*(2./d*atan(d/2.)));
378 
379  } else {
380  pq->Put(4, Zero3);
381  }
382 
383  // dot{q}
384  const Vec3& V0(pModalNode->GetVCurr());
385  const Vec3& W0(pModalNode->GetWCurr());
386  pqPrime->Put(1, RR.MulTV(V0));
387  pqPrime->Put(4, RR.MulTV(W0));
388 
389  // ddot{q}
390  pqPrime->Put(1, Zero3);
391  pqPrime->Put(4, Zero3);
392  }
393 
394  integer iModalIndex = pModalJoint->iGetModalIndex();
395 
396  for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
397  WorkVec.PutRowIndex(RigidF + iCnt, iModalIndex + NStModes + iCnt);
398  }
399 
400  integer iAeroIndex = iGetFirstIndex();
401  for (unsigned int iCnt = 1; iCnt <= NAeroStates + 2*NGust; iCnt++) {
402  WorkVec.PutRowIndex(RigidF + NStModes + iCnt, iAeroIndex + iCnt);
403  }
404 
405  /* Recupera i vettori {a} e {aP} e {aS} (deformate modali) */
406  for (unsigned int iCnt = 1; iCnt <= NStModes; iCnt++) {
407  pq->PutCoef(iCnt + RigidF, XCurr(iModalIndex + iCnt));
408  pqPrime->PutCoef(iCnt + RigidF, 0);
409  pqSec->PutCoef(iCnt + RigidF, 0);
410  }
411 
412  /* Recupera i vettori {xa} e {xaP} */
413  for (unsigned int iCnt = 1; iCnt <= NAeroStates; iCnt++) {
414  pxa->PutCoef(iCnt, XCurr(iAeroIndex + iCnt));
415  pxaPrime->PutCoef(iCnt, 0);
416  }
417 
418  for (unsigned int iCnt = 1; iCnt <= 2*NGust; iCnt++) {
419  pgs->PutCoef(iCnt, XCurr(iAeroIndex + NAeroStates + iCnt));
420  pgsPrime->PutCoef(iCnt, 0.);
421  }
422 
423  AssVec(WorkVec);
424 
425  return WorkVec;
426 }
427 
428 /* assemblaggio residuo */
429 void
430 AerodynamicModal::AssVec(SubVectorHandler& WorkVec)
431 {
432  DEBUGCOUTFNAME("AerodynamicModal::AssVec");
433 
434  /* Dati del nodo rigido */
435  const Vec3& X0(pModalNode->GetXCurr());
436  Vec3 V0(pModalNode->GetVCurr());
437  const Mat3x3& Rn(pModalNode->GetRCurr());
438  Mat3x3 RR(Rn*Ra);
439  Vec3 Vr(V0);
440 
441  Vr(1) = 0.0; // perturbazione di velocita', verificare
442 
443  doublereal rho, vs, p, T;
444  GetAirProps(X0, rho, vs, p, T); /* p, T are not used yet */
445 
446  Vec3 Vair(Zero3);
447  if (fGetAirVelocity(Vair, X0)) {
448  Vr -= Vair;
449  }
450  /* velocita' nel riferimento nodale aerodinamico */
451  V0 = RR.MulTV(Vr);
452  doublereal nV = std::abs(V0(1));
453  MyVectorHandler TmpA(NAeroStates);
454  TmpA.Reset();
455  MyVectorHandler Tmp(RigidF + NStModes);
456  Tmp.Reset();
457  MyVectorHandler TmpP(RigidF + NStModes);
458  TmpP.Reset();
459  MyVectorHandler TmpS(RigidF + NStModes);
460  TmpS.Reset();
461 
462  // FIXME: use matrix/vector product?
463  for (unsigned int i = 1; i <= NAeroStates; i++) {
464  for (unsigned int j = 1; j <= NStModes + RigidF; j++) {
465  TmpA.IncCoef(i, pB->operator()(i, j)*(pq->operator()(j)));
466  }
467  }
468 
469  pA->MatVecIncMul(TmpA, *pxa);
470 
471  /* doublereal CV = Chord/(2*nV); */
472  doublereal qd = 0.5*rho*nV*nV;
473  doublereal qd1 = 0.25*rho*nV*Chord; /* qd*CV */
474  doublereal qd2 = 0.125*rho*Chord*Chord; /* qd*CV*CV */
475 
476  for (unsigned int i = 1; i <= NAeroStates; i++) {
477  WorkVec.IncCoef(RigidF + NStModes + i, -pxaPrime->operator()(i) + (2*nV/Chord)*(TmpA(i)));
478  }
479 
480  for (unsigned int i = 1; i <= NStModes + RigidF; i++) {
481  for (unsigned int j = 1; j <= NStModes + RigidF; j++) {
482  Tmp.IncCoef(i, pD0->operator()(i, j)*(pq->operator()(j)));
483  TmpP.IncCoef(i, pD1->operator()(i, j)*(pqPrime->operator()(j)));
484  TmpS.IncCoef(i, pD2->operator()(i, j)*(pqSec->operator()(j)));
485  }
486  }
487 
488  for (unsigned int i = 1; i <= NStModes + RigidF; i++) {
489  for (unsigned int j = 1; j <= NAeroStates; j++) {
490  Tmp.IncCoef(i, pC->operator()(i, j)*(pxa->operator()(j)));
491  }
492  }
493 
494  if (RigidF) {
495  Vec3 F(Zero3);
496  Vec3 M(Zero3);
497 
498  for (unsigned int i = 1; i <= 3; i++) {
499  F.Put(i, qd*Tmp(i) + qd1*TmpP(i) + qd2*TmpS(i));
500  M.Put(i, qd*Tmp(i + 3) + qd1*TmpP(i + 3) + qd2*TmpS(i + 3));
501  }
502 
503  // std::cout << F << std::endl;
504  F = RR*F;
505  M = RR*M;
506 
507 #if 0
508  // tolti perche' altrimenti inizia un Dutch Roll
509  F(1) = 0.;
510  F(2) = 0.;
511  M(1) = 0.;
512  M(3) = 0.;
513 #endif
514 
515  WorkVec.Add(1, F);
516  WorkVec.Add(4, M);
517  }
518 
519  for (unsigned int i = 1 + RigidF; i <= NStModes + RigidF; i++) {
520  WorkVec.IncCoef(i, qd*Tmp(i) + qd1*TmpP(i) + qd2*TmpS(i));
521  }
522 
523  if (NGust) {
524 #if 0 /* unused */
525  doublereal Vyg = Vair(2)/nV;
526  doublereal Vzg = Vair(3)/nV;
527 #endif
528  for (unsigned int i = 1; i <= NAeroStates; i++) {
529  WorkVec.IncCoef(RigidF + NStModes + i,
530  (2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 1)*pgs->operator()(1)
531  + (2*nV/Chord)*pB->operator()(i, RigidF + NStModes + 2)*pgs->operator()(3));
532  }
533 
534  if (RigidF) {
535  Vec3 F(Zero3);
536  Vec3 M(Zero3);
537 
538  for (unsigned int i = 1; i <= 3; i++) {
539  F.Put(i, qd*pD0->operator()(i, RigidF + NStModes + 1)*pgs->operator()(1)
540  + qd*pD0->operator()(i, RigidF + NStModes + 2)*pgs->operator()(3)
541  + qd1*pD1->operator()(i, RigidF + NStModes + 1)*pgs->operator()(2)
542  + qd1*pD1->operator()(i, RigidF + NStModes + 2)*pgs->operator()(4)
543  + qd2*pD2->operator()(i, RigidF + NStModes + 1)*pgsPrime->operator()(2)
544  + qd2*pD2->operator()(i, RigidF + NStModes + 2)*pgsPrime->operator()(4));
545  M.Put(i, qd*pD0->operator()(i + 3, RigidF + NStModes + 1)*pgs->operator()(1)
546  + qd*pD0->operator()(i + 3, RigidF + NStModes + 2)*pgs->operator()(3)
547  + qd1*pD1->operator()(i + 3, RigidF + NStModes + 1)*pgs->operator()(2)
548  + qd1*pD1->operator()(i + 3, RigidF + NStModes + 2)*pgs->operator()(4)
549  + qd2*pD2->operator()(i + 3, RigidF + NStModes + 1)*pgsPrime->operator()(2)
550  + qd2*pD2->operator()(i + 3, RigidF + NStModes + 2)*pgsPrime->operator()(4));
551  }
552 
553  // std::cout << F << std::endl;
554  F = RR*F;
555  M = RR*M;
556  WorkVec.Add(1, F);
557  WorkVec.Add(4, M);
558  }
559 
560  for (unsigned int i = 1 + RigidF; i <= NStModes + RigidF; i++) {
561  WorkVec.IncCoef(i, qd*pD0->operator()(i, RigidF + NStModes + 1)*pgs->operator()(1)
562  + qd*pD0->operator()(i, RigidF + NStModes + 2)*pgs->operator()(3)
563  + qd1*pD1->operator()(i, RigidF + NStModes + 1)*pgs->operator()(2)
564  + qd1*pD1->operator()(i, RigidF + NStModes + 2)*pgs->operator()(4)
565  + qd2*pD2->operator()(i, RigidF + NStModes + 1)*pgsPrime->operator()(2)
566  + qd2*pD2->operator()(i, RigidF + NStModes + 2)*pgsPrime->operator()(4));
567  }
568 
569  for (unsigned int i = 0; i < NGust; i++) {
570  WorkVec.IncCoef(RigidF + i*2 + 1 + NStModes + NAeroStates,
571  -pgsPrime->operator()(1 + i*2) + pgs->operator()(2 + i*2));
572  if (nV != 0) {
573  WorkVec.IncCoef(RigidF + i*2 + 2 + NStModes + NAeroStates,
574  -pgsPrime->operator()(2 + i*2)
575  - gustVff*gustVff*pgs->operator()(1 + i*2)
576  - 2*gustXi*gustVff*pgs->operator()(2 + i*2)
577  + gustVff*gustVff*Vair(2 + i)/nV);
578 
579  } else {
580  WorkVec.IncCoef(RigidF + i*2 + 2 + NStModes + NAeroStates,
581  -pgsPrime->operator()(2 + i*2)
582  - gustVff*gustVff*pgs->operator()(1 + i*2)
583  - 2*gustXi*gustVff*pgs->operator()(2 + i*2));
584  }
585  }
586 
587  // std::cout << X0 << std::endl;
588  // std::cout << WorkVec(3) << std::endl;
589  }
590 }
591 
592 /* output; si assume che ogni tipo di elemento sappia, attraverso
593  * l'OutputHandler, dove scrivere il proprio output */
594 void
595 AerodynamicModal::Output(OutputHandler& OH) const
596 {
597  if (bToBeOutput()) {
598  OH.AeroModals() << std::setw(8) << GetLabel() << " ";
599  for (unsigned int iCnt = 1; iCnt <= NAeroStates; iCnt++) {
600  OH.AeroModals() << " " << pxa->operator()(iCnt);
601  }
602  for (unsigned int iCnt = 1; iCnt <= NAeroStates + 2*NGust; iCnt++) {
603  OH.AeroModals() << " " << pxaPrime->operator()(iCnt);
604  }
605  OH.AeroModals() << std::endl;
606  }
607 }
608 
609 /* AerodynamicModal - end */
610 
611 Elem *
612 ReadAerodynamicModal(DataManager* pDM,
613  MBDynParser& HP,
614  const DofOwner* pDO,
615  unsigned int uLabel)
616 {
617  DEBUGCOUTFNAME("ReadAerodynamicModal(" << uLabel << ")" << std::endl);
618 
619  /* formato dell'input:
620  *
621  * label,
622  * modal node,
623  * reference modal joint,
624  * rotation of the aerodynamic reference in respect of nodal reference,
625  * reference cord,
626  * number of aerodynamic states,
627  * {rigid, gust, gust filter cut-off frequency}
628  * file name containing state space model matrices;
629  */
630 
631  /* giunto modale collegato */
632  const Modal *pModalJoint = pDM->ReadElem<const Modal, const Joint, Elem::JOINT>(HP);
633  ReferenceFrame RF;
634  const StructNode *pModalNode = pModalJoint->pGetModalNode();
635  if (pModalNode) {
636  RF = ReferenceFrame(pModalNode);
637 
638  } else {
639  RF = ::AbsRefFrame;
640  }
641 
642  Mat3x3 Ra(HP.GetRotRel(RF));
643 
644  doublereal Chord = HP.GetReal();
645 
646  unsigned int AeroN = HP.GetInt();
647  /* numero modi e FEM */
648  unsigned int NModes = pModalJoint->uGetNModes();
649 
650  AerodynamicModal::RigidF_t rigidF = AerodynamicModal::NO_RIGID;
651  if (HP.IsKeyWord("rigid")) {
652  rigidF = AerodynamicModal::RIGID;
653  NModes += rigidF;
654  }
655 
656  /* Eventuale raffica */
657  unsigned int GustN = 0;
658  doublereal Vff = 0.;
659  if (HP.IsKeyWord("gust")) {
660  GustN = 2;
661  Vff = HP.GetReal();
662  }
663 
664  /* apre il file contenente le matrici A B C D0 D1 D2 */
665  const char *sFileData = HP.GetFileName();
666  std::ifstream fdat(sFileData);
667  DEBUGCOUT("Reading Aerodynamic State Space Matrices from file '"
668  << sFileData << '\'' << std::endl);
669  if (!fdat) {
670  int save_errno = errno;
671  silent_cerr(std::endl << "Unable to open file \"" << sFileData << "\" "
672  "(" << save_errno << ": " << strerror(save_errno) << ")" << std::endl);
673  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
674  }
675 
676  SpMapMatrixHandler* pAMat = 0;
677  FullMatrixHandler* pBMat = 0;
678  FullMatrixHandler* pCMat = 0;
679  FullMatrixHandler* pD0Mat = 0;
680  FullMatrixHandler* pD1Mat = 0;
681  FullMatrixHandler* pD2Mat = 0;
682  SAFENEWWITHCONSTRUCTOR(pAMat, SpMapMatrixHandler, SpMapMatrixHandler(AeroN, AeroN));
683  pAMat->Reset();
684  SAFENEWWITHCONSTRUCTOR(pBMat, FullMatrixHandler, FullMatrixHandler(AeroN, NModes + GustN));
685  SAFENEWWITHCONSTRUCTOR(pCMat, FullMatrixHandler, FullMatrixHandler(NModes, AeroN));
686  SAFENEWWITHCONSTRUCTOR(pD0Mat, FullMatrixHandler, FullMatrixHandler(NModes, NModes + GustN));
687  SAFENEWWITHCONSTRUCTOR(pD1Mat, FullMatrixHandler, FullMatrixHandler(NModes, NModes + GustN));
688  SAFENEWWITHCONSTRUCTOR(pD2Mat, FullMatrixHandler, FullMatrixHandler(NModes, NModes + GustN));
689 
690  doublereal d;
691  char str[BUFSIZ];
692 
693  /* parsing del file */
694  while (!fdat.eof()) {
695  fdat.getline(str, sizeof(str));
696 
697  /* legge il primo blocco (HEADER) */
698  if (!strncmp("*** MATRIX A", str, STRLENOF("*** MATRIX A"))) {
699  for (unsigned int iCnt = 1; iCnt <= AeroN; iCnt++) {
700  for (unsigned int jCnt = 1; jCnt <= AeroN; jCnt++) {
701  fdat >> d;
702  pAMat->PutCoef(iCnt, jCnt, d);
703  }
704  }
705 
706  /* legge il primo blocco (HEADER) */
707  } else if (!strncmp("*** MATRIX B", str, STRLENOF("*** MATRIX B"))) {
708  for (unsigned int iCnt = 1; iCnt <= AeroN; iCnt++) {
709  for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
710  fdat >> d;
711  pBMat->PutCoef(iCnt, jCnt, d);
712  }
713  }
714 
715  /* legge il primo blocco (HEADER) */
716  } else if (!strncmp("*** MATRIX C", str, STRLENOF("*** MATRIX C"))) {
717  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
718  for (unsigned int jCnt = 1; jCnt <= AeroN; jCnt++) {
719  fdat >> d;
720  pCMat->PutCoef(iCnt, jCnt, d);
721  }
722  }
723 
724  /* legge il primo blocco (HEADER) */
725  } else if (!strncmp("*** MATRIX D0", str, STRLENOF("*** MATRIX D0"))) {
726  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
727  for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
728  fdat >> d;
729  pD0Mat->PutCoef(iCnt, jCnt, d);
730  }
731  }
732 
733  /* legge il primo blocco (HEADER) */
734  } else if (!strncmp("*** MATRIX D1", str, STRLENOF("*** MATRIX D1"))) {
735  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
736  for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
737  fdat >> d;
738  pD1Mat->PutCoef(iCnt, jCnt, d);
739  }
740  }
741 
742  /* legge il primo blocco (HEADER) */
743  } else if (!strncmp("*** MATRIX D2", str, STRLENOF("*** MATRIX D2"))) {
744  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
745  for (unsigned int jCnt = 1; jCnt <= NModes + GustN; jCnt++) {
746  fdat >> d;
747  pD2Mat->PutCoef(iCnt, jCnt, d);
748  }
749  }
750  }
751  }
752  fdat.close();
753 
754  flag fOut = pDM->fReadOutput(HP, Elem::AEROMODAL);
755 
756  Elem* pEl = 0;
757  SAFENEWWITHCONSTRUCTOR(pEl,
758  AerodynamicModal,
759  AerodynamicModal(uLabel, pModalNode, pModalJoint,
760  Ra, pDO, Chord, NModes - rigidF,
761  AeroN, rigidF, GustN, Vff,
762  pAMat, pBMat, pCMat, pD0Mat, pD1Mat, pD2Mat, fOut));
763 
764  /* Se non c'e' il punto e virgola finale */
765  if (HP.IsArg()) {
766  silent_cerr("semicolon expected at line "
767  << HP.GetLineData() << std::endl);
768  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
769  }
770 
771  return pEl;
772 } /* End of DataManager::ReadAerodynamicModal() */
773 
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
AerodynamicModal::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: aeromodal.cc:278
MBDynParser::RF
RFType RF
Definition: mbpar.h:166
Modal::pGetModalNode
const ModalNode * pGetModalNode(void) const
Definition: modal.h:395
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
AerodynamicModal::pq
MyVectorHandler * pq
Definition: aeromodal.h:93
AerodynamicModal::pModalNode
const StructNode * pModalNode
Definition: aeromodal.h:75
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
Zero3
const Vec3 Zero3(0., 0., 0.)
CGR_Rot::Param
const Param_Manip Param
Definition: matvec3.cc:644
AerodynamicModal::AssVec
void AssVec(SubVectorHandler &WorkVec)
Definition: aeromodal.cc:430
AirPropOwner::GetAirProps
virtual bool GetAirProps(const Vec3 &X, doublereal &rho, doublereal &c, doublereal &p, doublereal &T) const
Definition: aerodyn.cc:760
AerodynamicModal::pB
FullMatrixHandler * pB
Definition: aeromodal.h:87
flag
long int flag
Definition: mbdyn.h:43
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
ReadAerodynamicModal
Elem * ReadAerodynamicModal(DataManager *pDM, MBDynParser &HP, const DofOwner *pDO, unsigned int uLabel)
Definition: aeromodal.cc:612
Vec3
Definition: matvec3.h:98
AerodynamicModal::WorkSpaceDim
void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: aeromodal.h:163
DEBUGCOUTFNAME
#define DEBUGCOUTFNAME(fname)
Definition: myassert.h:256
VectorHandler::ResizeReset
virtual void ResizeReset(integer)
Definition: vh.cc:55
AerodynamicModal::NStModes
unsigned int NStModes
Definition: aeromodal.h:82
HighParser::GetInt
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
AerodynamicModal::pC
FullMatrixHandler * pC
Definition: aeromodal.h:88
AerodynamicModal::pA
SpMapMatrixHandler * pA
Definition: aeromodal.h:86
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
AerodynamicModal::pD2
FullMatrixHandler * pD2
Definition: aeromodal.h:91
StructDispNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: strnode.cc:145
DataManager::ReadElem
Elem * ReadElem(MBDynParser &HP, Elem::Type type) const
Definition: dataman3.cc:2334
OutputHandler::AeroModals
std::ostream & AeroModals(void) const
Definition: output.h:527
AerodynamicModal::Ra
const Mat3x3 Ra
Definition: aeromodal.h:79
MatrixHandler::MatVecIncMul
virtual VectorHandler & MatVecIncMul(VectorHandler &out, const VectorHandler &in) const
Definition: mh.cc:350
Modal::uGetNModes
integer uGetNModes(void) const
Definition: modal.h:363
VectorHandler::IncCoef
virtual void IncCoef(integer iRow, const doublereal &dCoef)=0
AerodynamicModal::pgs
MyVectorHandler * pgs
Definition: aeromodal.h:100
IncludeParser::GetFileName
virtual const char * GetFileName(enum Delims Del=DEFAULTDELIM)
Definition: parsinc.cc:673
MyVectorHandler::Put
virtual void Put(integer iRow, const Vec3 &v)
Definition: vh.cc:699
AerodynamicModal::~AerodynamicModal
~AerodynamicModal(void)
Definition: aeromodal.cc:102
AerodynamicModal::pqSec
MyVectorHandler * pqSec
Definition: aeromodal.h:95
FullSubMatrixHandler::IncCoef
void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:683
AerodynamicModal::pModalJoint
const Modal * pModalJoint
Definition: aeromodal.h:76
AerodynamicModal::pxaPrime
MyVectorHandler * pxaPrime
Definition: aeromodal.h:98
MyVectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)
Definition: vh.h:261
RotManip::VecRot
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
AerodynamicModal::pD0
FullMatrixHandler * pD0
Definition: aeromodal.h:89
Mat3x3::MulTV
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
AerodynamicModal::RigidF
RigidF_t RigidF
Definition: aeromodal.h:113
Modal
Definition: modal.h:74
AerodynamicModal::pxa
MyVectorHandler * pxa
Definition: aeromodal.h:97
SubVectorHandler::PutRowIndex
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
SpMapMatrixHandler::Reset
void Reset(void)
Definition: spmapmh.cc:161
AbsRefFrame
const ReferenceFrame AbsRefFrame(0, Vec3(0., 0., 0), Mat3x3(1., 0., 0., 0., 1., 0., 0., 0., 1.), Vec3(0., 0., 0), Vec3(0., 0., 0), EULER_123)
AerodynamicModal::Restart
std::ostream & Restart(std::ostream &out) const
Definition: aeromodal.cc:142
FullMatrixHandler::PutCoef
virtual void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: fullmh.h:178
AerodynamicModal::gustVff
const doublereal gustVff
Definition: aeromodal.h:103
MBDynParser::pDM
DataManager * pDM
Definition: mbpar.h:252
AerodynamicModal::gustXi
const doublereal gustXi
Definition: aeromodal.h:104
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
Gust
Definition: gust.h:44
AerodynamicModal::InitialAssRes
SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: aeromodal.cc:350
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
AerodynamicModal::NAeroStates
unsigned int NAeroStates
Definition: aeromodal.h:83
StructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const =0
StructNode::GetWCurr
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
AerodynamicModal::Chord
doublereal Chord
Definition: aeromodal.h:81
StructNode::GetWPCurr
virtual const Vec3 & GetWPCurr(void) const
Definition: strnode.h:1042
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
MyVectorHandler::Reset
virtual void Reset(void)
Definition: vh.cc:459
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
VectorHandler::Add
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
AerodynamicModal::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &, const VectorHandler &)
Definition: aeromodal.cc:149
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
AerodynamicModal::Output
void Output(OutputHandler &OH) const
Definition: aeromodal.cc:595
AerodynamicModal::pqPrime
MyVectorHandler * pqPrime
Definition: aeromodal.h:94
SpMapMatrixHandler::PutCoef
void PutCoef(integer ix, integer iy, const doublereal &val)
Definition: spmapmh.h:187
Modal::iGetModalIndex
integer iGetModalIndex(void) const
Definition: modal.h:391
STRLENOF
#define STRLENOF(s)
Definition: mbdyn.h:166
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
AerodynamicModal::AerodynamicModal
AerodynamicModal(unsigned int uLabel, const StructNode *pN, const Modal *pMJ, const Mat3x3 &RaTmp, const DofOwner *pDO, doublereal Cd, const int NModal, const int NAero, RigidF_t rgF, const int Gust, const doublereal Vff, SpMapMatrixHandler *pAMat, FullMatrixHandler *pBMat, FullMatrixHandler *pCMat, FullMatrixHandler *pD0Mat, FullMatrixHandler *pD1Mat, FullMatrixHandler *pD2Mat, flag fout)
Definition: aeromodal.cc:43
Elem
Definition: elem.h:75
AerodynamicModal::pD1
FullMatrixHandler * pD1
Definition: aeromodal.h:90
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
AerodynamicModal::pgsPrime
MyVectorHandler * pgsPrime
Definition: aeromodal.h:101
AerodynamicModal::NGust
unsigned int NGust
Definition: aeromodal.h:84
StructDispNode::GetVCurr
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
AirPropOwner::fGetAirVelocity
virtual flag fGetAirVelocity(Vec3 &Velocity, const Vec3 &X) const
Definition: aerodyn.cc:725
Joint
Definition: joint.h:50
grad::atan
GradientExpression< UnaryExpr< FuncAtan, Expr > > atan(const GradientExpression< Expr > &u)
Definition: gradient.h:2985
StructDispNode::GetXPPCurr
virtual const Vec3 & GetXPPCurr(void) const
Definition: strnode.h:334
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Vec3::Put
void Put(unsigned short int iRow, const doublereal &dCoef)
Definition: matvec3.h:276
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056
MyVectorHandler::IncCoef
virtual void IncCoef(integer iRow, const doublereal &dCoef)
Definition: vh.h:278