MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
windprof.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/windprof.cc,v 1.13 2017/01/12 14:45:58 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /*
33  * Wind profiles
34  */
35 
36 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
37 
38 #include "dataman.h"
39 #include "mbpar.h"
40 
41 #include "windprof.h"
42 
43 /* WindProfile - begin */
44 
46  const Vec3& X0,
47  const Mat3x3& R0)
48 : X0(X0), R0(R0)
49 {
50  NO_OP;
51 }
52 
54 {
55  NO_OP;
56 }
57 
58 /* WindProfile - end */
59 
60 /* ScalarFuncWindProfile - begin */
61 
63  const Vec3& X0,
64  const Mat3x3& R0,
65  const BasicScalarFunction *sf)
66 : WindProfile(X0, R0),
67 sf(sf)
68 {
69  ASSERT(sf != 0);
70 }
71 
73 {
74  NO_OP;
75 }
76 
77 bool
79 {
80  // dZ is the component of (X - X0) along axis 3 of R0
81  doublereal dZ = R0.GetVec(3)*(X - X0);
82 
83  // V is projected along axis 1 of R0
84  V = R0.GetVec(1)*(*sf)(dZ);
85 
86  return true;
87 }
88 
89 std::ostream&
90 ScalarFuncWindProfile::Restart(std::ostream& out) const
91 {
92  return out << "scalar function, "
93  << ", reference position, " << X0
94  << ", reference orientation, " << R0
95  << " # not implemented yet";
96 }
97 
99 {
100  NO_OP;
101 }
102 
103 Gust *
105 {
106  const BasicScalarFunction *sf = 0;
107  Vec3 X0(Zero3);
108  bool bGotX0 = false;
109  Mat3x3 R0(Eye3);
110  bool bGotR0 = false;
111 
112  while (HP.IsArg()) {
113  if (HP.IsKeyWord("reference" "position")) {
114  if (bGotX0) {
115  silent_cerr("ScalarFuncWindProfile: "
116  "reference position provided twice "
117  "at line " << HP.GetLineData()
118  << std::endl);
120  }
121 
122  X0 = HP.GetVecAbs(::AbsRefFrame);
123 
124  } else if (HP.IsKeyWord("reference" "orientation")) {
125  if (bGotR0) {
126  silent_cerr("ScalarFuncWindProfile: "
127  "reference orientation provided twice "
128  "at line " << HP.GetLineData()
129  << std::endl);
131  }
132 
133  R0 = HP.GetRotAbs(::AbsRefFrame);
134 
135  } else {
136  break;
137  }
138  }
139 
140  if (!HP.IsArg()) {
141  silent_cerr("ScalarFuncWindProfile: "
142  "function expected "
143  "at line " << HP.GetLineData() << std::endl);
145  }
146 
147  sf = ParseScalarFunction(HP, const_cast<DataManager *const>(pDM));
148  if (sf == 0) {
149  silent_cerr("ScalarFuncWindProfile: "
150  "unable to parse scalar function "
151  "at line " << HP.GetLineData() << std::endl);
153  }
154 
155  Gust *pG = 0;
157  ScalarFuncWindProfile(X0, R0, sf));
158 
159  return pG;
160 }
161 
162 /* ScalarFuncWindProfile - end */
163 
164 /* PowerLawWindProfile - begin */
165 
167  const Vec3& X0,
168  const Mat3x3& R0,
169  const doublereal dZRef,
170  const DriveCaller *pVRef,
171  const doublereal dPower)
172 : WindProfile(X0, R0),
173 dZRef(dZRef),
174 VRef(pVRef),
175 dPower(dPower)
176 {
177  ASSERT(dZRef > 0.);
178  ASSERT(pVRef != 0);
179  // NOTE: should be < 1?
180  ASSERT(dPower > 0.);
181 }
182 
184 {
185  NO_OP;
186 }
187 
188 bool
190 {
191  // dZ is the component of (X - X0) along axis 3 of R0
192  doublereal dZ = R0.GetVec(3)*(X - X0);
193  if (dZ <= 0.) {
194  return false;
195  }
196 
197  // V is projected along axis 1 of R0
198  V = R0.GetVec(1)*(VRef.dGet()*std::pow(dZ/dZRef, dPower));
199 
200  return true;
201 }
202 
203 std::ostream&
204 PowerLawWindProfile::Restart(std::ostream& out) const
205 {
206  out << "power law"
207  << ", reference position, " << X0
208  << ", reference orientation, " << R0
209  << ", reference elevation, " << dZRef
210  << ", reference velocity, ", VRef.pGetDriveCaller()->Restart(out)
211  << ", exponent, " << dPower;
212  return out;
213 }
214 
216 {
217  NO_OP;
218 }
219 
220 Gust *
222 {
223  doublereal dZRef = -1.;
224  DriveCaller *pVRef = 0;
225  doublereal dPower = -1.;
226  Vec3 X0(Zero3);
227  bool bGotX0 = false;
228  Mat3x3 R0(Eye3);
229  bool bGotR0 = false;
230 
231  while (HP.IsArg()) {
232  if (HP.IsKeyWord("reference" "elevation")) {
233  if (dZRef != -1.) {
234  silent_cerr("PowerLawWindProfile: "
235  "reference elevation provided twice "
236  "at line " << HP.GetLineData()
237  << std::endl);
239  }
240 
241  dZRef = HP.GetReal();
242 
243  } else if (HP.IsKeyWord("reference" "velocity")) {
244  if (pVRef != 0) {
245  silent_cerr("PowerLawWindProfile: "
246  "reference velocity provided twice "
247  "at line " << HP.GetLineData()
248  << std::endl);
250  }
251 
252  pVRef = HP.GetDriveCaller();
253 
254  } else if (HP.IsKeyWord("exponent")) {
255  if (dPower != -1.) {
256  silent_cerr("PowerLawWindProfile: "
257  "exponent provided twice "
258  "at line " << HP.GetLineData()
259  << std::endl);
261  }
262 
263  dPower = HP.GetReal();
264 
265  } else if (HP.IsKeyWord("reference" "position")) {
266  if (bGotX0) {
267  silent_cerr("PowerLawWindProfile: "
268  "reference position provided twice "
269  "at line " << HP.GetLineData()
270  << std::endl);
272  }
273 
274  X0 = HP.GetVecAbs(::AbsRefFrame);
275 
276  } else if (HP.IsKeyWord("reference" "orientation")) {
277  if (bGotR0) {
278  silent_cerr("PowerLawWindProfile: "
279  "reference orientation provided twice "
280  "at line " << HP.GetLineData()
281  << std::endl);
283  }
284 
285  R0 = HP.GetRotAbs(::AbsRefFrame);
286 
287  } else {
288  break;
289  }
290  }
291 
292  if (dZRef == -1.) {
293  silent_cerr("PowerLawWindProfile: "
294  "reference elevation missing "
295  "at line " << HP.GetLineData()
296  << std::endl);
298  }
299 
300  if (pVRef == 0) {
301  silent_cerr("PowerLawWindProfile: "
302  "reference velocity missing "
303  "at line " << HP.GetLineData()
304  << std::endl);
306  }
307 
308  if (dPower == -1.) {
309  silent_cerr("PowerLawWindProfile: "
310  "exponent missing "
311  "at line " << HP.GetLineData()
312  << std::endl);
314  }
315 
316  Gust *pG = 0;
318  PowerLawWindProfile(X0, R0, dZRef, pVRef, dPower));
319 
320  return pG;
321 }
322 
323 /* PowerLawWindProfile - end */
324 
325 /* LogarithmicWindProfile - begin */
326 
328  const Vec3& X0,
329  const Mat3x3& R0,
330  const doublereal dZRef,
331  const DriveCaller *pVRef,
332  const doublereal dSurfaceRoughnessLength)
333 : WindProfile(X0, R0),
334 dZRef(dZRef),
335 VRef(pVRef),
336 dSurfaceRoughnessLength(dSurfaceRoughnessLength),
337 logZRefZ0(std::log(dZRef/dSurfaceRoughnessLength))
338 {
339  ASSERT(dZRef > 0.);
340  ASSERT(pVRef != 0);
341  ASSERT(dSurfaceRoughnessLength > 0.);
342 }
343 
345 {
346  NO_OP;
347 }
348 
349 bool
351 {
352  // dZ is the component of (X - X0) along axis 3 of R0
353  doublereal dZ = R0.GetVec(3)*(X - X0);
354  if (dZ <= 0.) {
355  return false;
356  }
357 
358  // V is projected along axis 1 of R0
359  V = R0.GetVec(1)*(VRef.dGet()*(std::log(dZ/dSurfaceRoughnessLength) - 0.)/(logZRefZ0 - 0.));
360 
361  return true;
362 }
363 
364 std::ostream&
365 LogarithmicWindProfile::Restart(std::ostream& out) const
366 {
367  out << "logarithmic"
368  << ", reference position, " << X0
369  << ", reference orientation, " << R0
370  << ", reference elevation, " << dZRef
371  << ", reference velocity, ", VRef.pGetDriveCaller()->Restart(out)
372  << ", surface roughness length, " << dSurfaceRoughnessLength;
373  return out;
374 }
375 
377 {
378  NO_OP;
379 }
380 
381 Gust *
383 {
384  doublereal dZRef = -1.;
385  DriveCaller *pVRef = 0;
386  doublereal dSurfaceRoughnessLength = -1.;
387  Vec3 X0(Zero3);
388  bool bGotX0 = false;
389  Mat3x3 R0(Eye3);
390  bool bGotR0 = false;
391 
392  while (HP.IsArg()) {
393  if (HP.IsKeyWord("reference" "elevation")) {
394  if (dZRef != -1.) {
395  silent_cerr("LogarithmicWindProfile: "
396  "reference elevation provided twice "
397  "at line " << HP.GetLineData()
398  << std::endl);
400  }
401 
402  dZRef = HP.GetReal();
403 
404  } else if (HP.IsKeyWord("reference" "velocity")) {
405  if (pVRef != 0) {
406  silent_cerr("LogarithmicWindProfile: "
407  "reference velocity provided twice "
408  "at line " << HP.GetLineData()
409  << std::endl);
411  }
412 
413  pVRef = HP.GetDriveCaller();
414 
415  } else if (HP.IsKeyWord("surface" "roughness" "length")) {
416  if (dSurfaceRoughnessLength != -1.) {
417  silent_cerr("LogarithmicWindProfile: "
418  "surface roughness length provided twice "
419  "at line " << HP.GetLineData()
420  << std::endl);
422  }
423 
424  dSurfaceRoughnessLength = HP.GetReal();
425 
426  } else if (HP.IsKeyWord("reference" "position")) {
427  if (bGotX0) {
428  silent_cerr("PowerLawWindProfile: "
429  "reference position provided twice "
430  "at line " << HP.GetLineData()
431  << std::endl);
433  }
434 
435  X0 = HP.GetVecAbs(::AbsRefFrame);
436 
437  } else if (HP.IsKeyWord("reference" "orientation")) {
438  if (bGotR0) {
439  silent_cerr("PowerLawWindProfile: "
440  "reference orientation provided twice "
441  "at line " << HP.GetLineData()
442  << std::endl);
444  }
445 
446  R0 = HP.GetRotAbs(::AbsRefFrame);
447 
448  } else {
449  break;
450  }
451  }
452 
453  if (dZRef == -1.) {
454  silent_cerr("LogarithmicWindProfile: "
455  "reference elevation missing "
456  "at line " << HP.GetLineData()
457  << std::endl);
459  }
460 
461  if (pVRef == 0) {
462  silent_cerr("LogarithmicWindProfile: "
463  "reference velocity missing "
464  "at line " << HP.GetLineData()
465  << std::endl);
467  }
468 
469  if (dSurfaceRoughnessLength == -1.) {
470  silent_cerr("LogarithmicWindProfile: "
471  "surface roughness length missing "
472  "at line " << HP.GetLineData()
473  << std::endl);
475  }
476 
477  Gust *pG = 0;
479  LogarithmicWindProfile(X0, R0,
480  dZRef, pVRef, dSurfaceRoughnessLength));
481 
482  return pG;
483 }
484 
485 /* LogarithmicWindProfile - end */
486 
virtual ~PowerLawWindProfile(void)
Definition: windprof.cc:183
const Vec3 Zero3(0., 0., 0.)
virtual Gust * Read(const DataManager *pDM, MBDynParser &HP)
Definition: windprof.cc:382
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
DriveOwner VRef
Definition: windprof.h:87
Mat3x3 GetRotAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1857
virtual Gust * Read(const DataManager *pDM, MBDynParser &HP)
Definition: windprof.cc:221
const doublereal dZRef
Definition: windprof.h:86
virtual bool GetVelocity(const Vec3 &X, Vec3 &V) const
Definition: windprof.cc:350
virtual std::ostream & Restart(std::ostream &out) const
Definition: windprof.cc:90
virtual Gust * Read(const DataManager *pDM, MBDynParser &HP)
Definition: windprof.cc:104
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
#define NO_OP
Definition: myassert.h:74
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
virtual std::ostream & Restart(std::ostream &out) const =0
const doublereal logZRefZ0
Definition: windprof.h:116
const doublereal dPower
Definition: windprof.h:88
virtual ~LogarithmicGR(void)
Definition: windprof.cc:376
LogarithmicWindProfile(const Vec3 &X0, const Mat3x3 &R0, const doublereal dZRef, const DriveCaller *pVRef, const doublereal dSurfaceRoughnessLength)
Definition: windprof.cc:327
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)
GradientExpression< UnaryExpr< FuncLog, Expr > > log(const GradientExpression< Expr > &u)
Definition: gradient.h:2976
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
ScalarFuncWindProfile(const Vec3 &X0, const Mat3x3 &R0, const BasicScalarFunction *sf)
Definition: windprof.cc:62
Definition: gust.h:44
virtual ~ScalarFuncGR(void)
Definition: windprof.cc:98
virtual ~PowerLawGR(void)
Definition: windprof.cc:215
WindProfile(const Vec3 &X0, const Mat3x3 &R0)
Definition: windprof.cc:45
virtual ~ScalarFuncWindProfile(void)
Definition: windprof.cc:72
#define ASSERT(expression)
Definition: colamd.c:977
const Mat3x3 R0
Definition: windprof.h:50
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
virtual std::ostream & Restart(std::ostream &out) const
Definition: windprof.cc:204
Vec3 GetVecAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1641
virtual ~WindProfile(void)
Definition: windprof.cc:53
virtual bool IsArg(void)
Definition: parser.cc:807
virtual bool GetVelocity(const Vec3 &X, Vec3 &V) const
Definition: windprof.cc:189
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
virtual ~LogarithmicWindProfile(void)
Definition: windprof.cc:344
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
virtual std::ostream & Restart(std::ostream &out) const
Definition: windprof.cc:365
virtual bool GetVelocity(const Vec3 &X, Vec3 &V) const
Definition: windprof.cc:78
double doublereal
Definition: colamd.c:52
const doublereal dSurfaceRoughnessLength
Definition: windprof.h:114
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
const BasicScalarFunction *const ParseScalarFunction(MBDynParser &HP, DataManager *const pDM)
const doublereal dZRef
Definition: windprof.h:112
PowerLawWindProfile(const Vec3 &X0, const Mat3x3 &R0, const doublereal dZRef, const DriveCaller *pVRef, const doublereal dPower)
Definition: windprof.cc:166
const Vec3 X0
Definition: windprof.h:49
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056