Nektar++
Static Public Member Functions | Static Public Attributes | Protected Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
Nektar::ForcingMovingBody Class Reference

#include <ForcingMovingBody.h>

Inheritance diagram for Nektar::ForcingMovingBody:
Inheritance graph
[legend]
Collaboration diagram for Nektar::ForcingMovingBody:
Collaboration graph
[legend]

Static Public Member Functions

static SolverUtils::ForcingSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
 Creates an instance of this class. More...
 
- Static Public Member Functions inherited from Nektar::SolverUtils::Forcing
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
 

Static Public Attributes

static std::string className
 Name of the class. More...
 

Protected Member Functions

virtual void v_InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
 
virtual void v_Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
 
- Protected Member Functions inherited from Nektar::SolverUtils::Forcing
SOLVER_UTILS_EXPORT Forcing (const LibUtilities::SessionReaderSharedPtr &)
 Constructor. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (Array< OneD, MultiRegions::ExpListSharedPtr > pFields, LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
 
SOLVER_UTILS_EXPORT void EvaluateTimeFunction (LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
 

Private Member Functions

 ForcingMovingBody (const LibUtilities::SessionReaderSharedPtr &pSession)
 
void CheckIsFromFile (const TiXmlElement *pForce)
 
void InitialiseCableModel (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 
void InitialiseFilter (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
 
void UpdateMotion (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
 
void TensionedCableModel (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
 
void EvaluateStructDynModel (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, NekDouble time)
 
void CalculateForcing (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
void MappingBndConditions (const Array< OneD, MultiRegions::ExpListSharedPtr > &pfields, const Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
 
void EvaluateAccelaration (const Array< OneD, NekDouble > &input, Array< OneD, NekDouble > &output, int npoints)
 
void SetDynEqCoeffMatrix (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 
void RollOver (Array< OneD, Array< OneD, NekDouble > > &input)
 

Private Attributes

int m_movingBodyCalls
 number of times the movbody have been called More...
 
int m_np
 number of planes per processors More...
 
int m_vdim
 vibration dimension More...
 
NekDouble m_structrho
 mass of the cable per unit length More...
 
NekDouble m_structdamp
 damping ratio of the cable More...
 
NekDouble m_lhom
 length ratio of the cable More...
 
NekDouble m_kinvis
 fluid viscous More...
 
NekDouble m_timestep
 time step More...
 
LibUtilities::NektarFFTSharedPtr m_FFT
 
FilterMovingBodySharedPtr m_MovBodyfilter
 
Array< OneD, NekDouble > m_Aeroforces
 storage for the cable's force(x,y) variables More...
 
Array< OneD, NekDouble > m_MotionVars
 storage for the cable's motion(x,y) variables More...
 
Array< OneD, Array< OneD, NekDouble > > m_W
 srorage for the velocity in z-direction More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
 fictitious velocity storage More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
 fictitious acceleration storage More...
 
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
 matrices in Newmart-beta method More...
 
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
 matrices in Newmart-beta method More...
 
Array< OneD, std::string > m_funcName
 [0] is displacements, [1] is velocities, [2] is accelerations More...
 
Array< OneD, std::string > m_motion
 motion direction: [0] is 'x' and [1] is 'y' More...
 
Array< OneD, bool > m_IsFromFile
 do determine if the the body motion come from an extern file More...
 
Array< OneD, Array< OneD, NekDouble > > m_zta
 Store the derivatives of motion variables in x-direction. More...
 
Array< OneD, Array< OneD, NekDouble > > m_eta
 Store the derivatives of motion variables in y-direction. More...
 
Array< OneD, Array< OneD, NekDouble > > m_forcing
 

Friends

class MemoryManager< ForcingMovingBody >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Forcing
virtual SOLVER_UTILS_EXPORT ~Forcing ()
 
SOLVER_UTILS_EXPORT void InitObject (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
 Initialise the forcing object. More...
 
SOLVER_UTILS_EXPORT void Apply (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
 Apply the forcing. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::Forcing
LibUtilities::SessionReaderSharedPtr m_session
 Session reader. More...
 
Array< OneD, Array< OneD, NekDouble > > m_Forcing
 Evaluated forcing function. More...
 
int m_NumVariable
 Number of variables. More...
 

Detailed Description

Definition at line 49 of file ForcingMovingBody.h.

Constructor & Destructor Documentation

Nektar::ForcingMovingBody::ForcingMovingBody ( const LibUtilities::SessionReaderSharedPtr &  pSession)
private

Definition at line 47 of file ForcingMovingBody.cpp.

49  : Forcing(pSession)
50 {
51 }
Nektar::SolverUtils::Forcing::Forcing
SOLVER_UTILS_EXPORT Forcing(const LibUtilities::SessionReaderSharedPtr &)
Constructor.
Definition: Forcing.cpp:51

Member Function Documentation

void Nektar::ForcingMovingBody::CalculateForcing ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields)
private

Definition at line 245 of file ForcingMovingBody.cpp.

References Nektar::MultiRegions::DirCartesianMap, EvaluateAccelaration(), m_eta, m_forcing, m_kinvis, m_zta, Vmath::Smul(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Zero().

Referenced by v_Apply().

247 {
248 
249  int nPointsTot = fields[0]->GetNpoints();
250  Array<OneD, NekDouble> U,V,W;
251  Array<OneD, NekDouble> Wt,Wx,Wy,Wz,Wxx,Wxy,Wxz,Wyy,Wyz,Wzz;
252  Array<OneD, NekDouble> Px,Py,Pz;
253  Array<OneD, NekDouble> tmp,tmp1,tmp2,tmp3;
254  Array<OneD, NekDouble> Fx,Fy,Fz;
255  U = Array<OneD, NekDouble> (nPointsTot);
256  V = Array<OneD, NekDouble> (nPointsTot);
257  W = Array<OneD, NekDouble> (nPointsTot);
258  Wt = Array<OneD, NekDouble> (nPointsTot);
259  Wx = Array<OneD, NekDouble> (nPointsTot);
260  Wy = Array<OneD, NekDouble> (nPointsTot);
261  Wz = Array<OneD, NekDouble> (nPointsTot);
262  Wxx = Array<OneD, NekDouble> (nPointsTot);
263  Wxy = Array<OneD, NekDouble> (nPointsTot);
264  Wyy = Array<OneD, NekDouble> (nPointsTot);
265  Wxz = Array<OneD, NekDouble> (nPointsTot);
266  Wyz = Array<OneD, NekDouble> (nPointsTot);
267  Wzz = Array<OneD, NekDouble> (nPointsTot);
268  Px = Array<OneD, NekDouble> (nPointsTot);
269  Py = Array<OneD, NekDouble> (nPointsTot);
270  Pz = Array<OneD, NekDouble> (nPointsTot);
271  Fx = Array<OneD, NekDouble> (nPointsTot,0.0);
272  Fy = Array<OneD, NekDouble> (nPointsTot,0.0);
273  Fz = Array<OneD, NekDouble> (nPointsTot,0.0);
274  tmp = Array<OneD, NekDouble> (nPointsTot,0.0);
275  tmp1 = Array<OneD, NekDouble> (nPointsTot);
276  tmp2 = Array<OneD, NekDouble> (nPointsTot);
277  tmp3 = Array<OneD, NekDouble> (nPointsTot);
278  fields[0]->HomogeneousBwdTrans(fields[0]->GetPhys(),U);
279  fields[0]->HomogeneousBwdTrans(fields[1]->GetPhys(),V);
280  fields[0]->HomogeneousBwdTrans(fields[2]->GetPhys(),W);
281  fields[0]->HomogeneousBwdTrans(fields[3]->GetPhys(),tmp1); // pressure
282 
283  //-------------------------------------------------------------------------
284  // Setting the pressure derivatives
285  //-------------------------------------------------------------------------
286  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],tmp1,Px); // Px
287  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],tmp1,Py); // Py
288 
289  //-------------------------------------------------------------------------
290  // Setting the z-component velocity derivatives
291  //-------------------------------------------------------------------------
292  EvaluateAccelaration(W,Wt,nPointsTot); //Wt
293  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], W, Wx); // Wx
294  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], W, Wy); // Wy
295  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
296  fields[2]->GetPhys(),tmp1); // Wz
297  fields[0]->HomogeneousBwdTrans(tmp1, Wz); // Wz in physical space
298  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], Wx, Wxx); // Wxx
299  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], Wx, Wxy); // Wxy
300  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], Wy, Wyy); // Wyy
301  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], Wz, Wxz); // Wxz
302  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], Wz, Wyz); // Wyz
303  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2], tmp1, tmp2); // Wzz
304  fields[0]->HomogeneousBwdTrans(tmp2, Wzz); // Wzz in physical space
305  //-------------------------------------------------------------------------
306  // Setting the z-component of the accelaration term:
307  // tmp = (Wt + U * Wx + V * Wy + W * Wz)
308  //-------------------------------------------------------------------------
309  Vmath::Vadd(nPointsTot, tmp, 1, Wt, 1, tmp, 1);
310  Vmath::Vmul(nPointsTot, U, 1, Wx, 1, tmp1, 1);
311  Vmath::Vadd(nPointsTot, tmp, 1, tmp1, 1, tmp, 1);
312  Vmath::Vmul(nPointsTot, V, 1, Wy, 1, tmp1, 1);
313  Vmath::Vadd(nPointsTot, tmp, 1, tmp1, 1, tmp, 1);
314  Vmath::Vmul(nPointsTot, W, 1, Wz, 1, tmp1, 1);
315  Vmath::Vadd(nPointsTot, tmp, 1, tmp1, 1, tmp, 1);
316 
317  //-------------------------------------------------------------------------
318  // x-component of the m_forcing - accelaration component
319  //-------------------------------------------------------------------------
320  Vmath::Vsub(nPointsTot, Fx, 1, m_zta[2], 1, Fx, 1);
321  Vmath::Vmul(nPointsTot, m_zta[3], 1, tmp,1, tmp1, 1);
322  Vmath::Vsub(nPointsTot, Fx, 1, tmp1, 1, Fx, 1);
323 
324  Vmath::Vmul(nPointsTot, m_zta[6], 1, W, 1, tmp1, 1);
325  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
326  Vmath::Vsub(nPointsTot, Fx, 1, tmp2, 1, Fx, 1);
327  // W^2 - we reuse it later
328  Vmath::Vmul(nPointsTot, W, 1, W, 1, tmp3, 1);
329  Vmath::Vmul(nPointsTot, m_zta[4], 1, tmp3,1, tmp2, 1);
330  Vmath::Vsub(nPointsTot, Fx, 1, tmp2, 1, Fx, 1);
331 
332  //-------------------------------------------------------------------------
333  // y-component of the m_forcing - accelaration component
334  //-------------------------------------------------------------------------
335  Vmath::Vmul(nPointsTot, m_eta[4], 1, tmp3, 1, tmp2, 1); // reusing W^2
336  Vmath::Vsub(nPointsTot, Fy, 1, tmp2, 1, Fy, 1);
337 
338  Vmath::Vmul(nPointsTot, m_eta[3], 1, tmp, 1, tmp1, 1);
339  Vmath::Vsub(nPointsTot, Fy, 1, tmp1, 1, Fy, 1);
340 
341  Vmath::Vsub(nPointsTot, Fy, 1, m_eta[2], 1, Fy, 1);
342  Vmath::Vmul(nPointsTot, m_eta[6], 1, W, 1, tmp1, 1);
343  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
344  Vmath::Vsub(nPointsTot, Fy, 1, tmp2, 1, Fy, 1);
345 
346  //-------------------------------------------------------------------------
347  // z-component of the m_forcing - accelaration component
348  //-------------------------------------------------------------------------
349  Vmath::Vmul(nPointsTot, m_zta[3], 1, Px, 1, tmp1, 1);
350  Vmath::Vmul(nPointsTot, m_eta[3], 1, Py, 1, tmp2, 1);
351  Vmath::Vadd(nPointsTot, Fz, 1, tmp1, 1, Fz, 1);
352  Vmath::Vadd(nPointsTot, Fz, 1, tmp2, 1, Fz, 1);
353 
354  //-------------------------------------------------------------------------
355  // Note: Now we use Px,Py,Pz to store the viscous component of the m_forcing
356  // before we multiply them by m_kinvis = 1/Re. Since we build up on them we
357  // need to set the entries to zero.
358  //-------------------------------------------------------------------------
359  Vmath::Zero(nPointsTot,Px,1);
360  Vmath::Zero(nPointsTot,Py,1);
361  Vmath::Zero(nPointsTot,Pz,1);
362 
363  //-------------------------------------------------------------------------
364  // x-component of the m_forcing - viscous component1: (U_z^'z^' - U_zz + ZETA_tz^'z^')
365  //-------------------------------------------------------------------------
366  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
367  fields[0]->GetPhys(), tmp1); // Uz
368  fields[0]->HomogeneousBwdTrans(tmp1, tmp3); // Uz in physical space
369  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp3, tmp1); // Uzx
370  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], tmp3, tmp2); // Uzy
371 
372  Vmath::Vmul(nPointsTot, m_zta[3], 1, tmp1, 1, tmp1, 1);
373  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp1, 1);
374  Vmath::Vmul(nPointsTot, m_eta[3], 1, tmp2, 1, tmp2, 1);
375  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp2, 1);
376  Vmath::Vsub(nPointsTot, Px, 1, tmp1, 1, Px, 1);
377  Vmath::Vsub(nPointsTot, Px, 1, tmp2, 1, Px, 1);
378 
379  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], U, tmp1); // Ux
380  Vmath::Vmul(nPointsTot, m_zta[4], 1, tmp1, 1, tmp2, 1);
381  Vmath::Vsub(nPointsTot, Px, 1, tmp2, 1, Px, 1);
382  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp1, tmp2); // Uxx
383  Vmath::Vmul(nPointsTot, m_zta[9], 1, tmp2, 1, tmp3, 1);
384  Vmath::Vadd(nPointsTot, Px, 1, tmp3, 1, Px, 1);
385  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], tmp1, tmp2); // Uxy
386  Vmath::Vmul(nPointsTot, m_zta[8], 1, tmp2, 1, tmp3, 1);
387  Vmath::Smul(nPointsTot, 2.0, tmp3, 1, tmp3, 1);
388  Vmath::Vadd(nPointsTot, Px, 1, tmp3, 1, Px, 1);
389 
390  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], U, tmp1); // Uy
391  Vmath::Vmul(nPointsTot, m_eta[4], 1, tmp1, 1, tmp2, 1);
392  Vmath::Vsub(nPointsTot, Px, 1, tmp2, 1, Px, 1);
393  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], tmp1, tmp2); // Uyy
394  Vmath::Vmul(nPointsTot, m_eta[9], 1, tmp2, 1, tmp3, 1);
395  Vmath::Vadd(nPointsTot, Px, 1, tmp3, 1, Px, 1);
396  Vmath::Vadd(nPointsTot, Px, 1, m_zta[7], 1, Px, 1);
397 
398  //-------------------------------------------------------------------------
399  // x-component of the m_forcing - viscous component2:
400  // ((ZETA_z * W)_z^'z^' + ZETA_z * (W_xx + W_yy))
401  //-------------------------------------------------------------------------
402  Vmath::Vmul(nPointsTot, m_zta[5], 1, W, 1, tmp1, 1);
403  Vmath::Vadd(nPointsTot, Px, 1, tmp1, 1, Px, 1);
404 
405  Vmath::Vmul(nPointsTot, m_zta[3], 1, Wx, 1, tmp1, 1);
406  Vmath::Vmul(nPointsTot, m_zta[4], 1, tmp1, 1, tmp2, 1);
407  Vmath::Smul(nPointsTot, 3.0, tmp2, 1, tmp3, 1);
408  Vmath::Vsub(nPointsTot, Px, 1, tmp3, 1, Px, 1);
409 
410  Vmath::Vmul(nPointsTot, m_eta[3], 1, Wy, 1, tmp1, 1);
411  Vmath::Vmul(nPointsTot, m_zta[4], 1, tmp1, 1, tmp2, 1);
412  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp3, 1);
413  Vmath::Vsub(nPointsTot, Px, 1, tmp3, 1, Px, 1);
414 
415  Vmath::Vmul(nPointsTot, m_zta[3], 1, Wy, 1, tmp1, 1);
416  Vmath::Vmul(nPointsTot, m_eta[4], 1, tmp1, 1, tmp2, 1);
417  Vmath::Vsub(nPointsTot, Px, 1, tmp2, 1, Px, 1);
418 
419  Vmath::Vmul(nPointsTot, m_zta[4], 1, Wz, 1, tmp1, 1);
420  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
421  Vmath::Vadd(nPointsTot, Px, 1, tmp2, 1, Px, 1);
422 
423  Vmath::Vmul(nPointsTot, m_zta[9], 1, Wxy, 1, tmp1, 1);
424  Vmath::Vmul(nPointsTot, m_eta[3], 1, tmp1, 1, tmp2, 1);
425  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp3, 1);
426  Vmath::Vadd(nPointsTot, Px, 1, tmp3, 1, Px, 1);
427 
428  Vmath::Vmul(nPointsTot, m_zta[9], 1, Wxz, 1, tmp1, 1);
429  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
430  Vmath::Vsub(nPointsTot, Px, 1, tmp2, 1, Px, 1);
431 
432  Vmath::Vmul(nPointsTot, m_zta[8], 1, Wyz, 1, tmp1, 1);
433  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
434  Vmath::Vsub(nPointsTot, Px, 1, tmp2, 1, Px, 1);
435 
436  Vmath::Vmul(nPointsTot, m_zta[9], 1, m_zta[3], 1, tmp1, 1);
437  Vmath::Vmul(nPointsTot, tmp1,1, Wxx, 1, tmp2, 1);
438  Vmath::Vadd(nPointsTot, Px, 1, tmp2, 1, Px, 1);
439 
440  Vmath::Vmul(nPointsTot, m_zta[3], 1, Wyy, 1, tmp1, 1);
441  Vmath::Vmul(nPointsTot, m_eta[9], 1, tmp1, 1, tmp2, 1);
442  Vmath::Vadd(nPointsTot, Px, 1, tmp2, 1, Px, 1);
443 
444  Vmath::Vadd(nPointsTot, Wxx, 1, Wyy, 1, tmp1, 1);
445  Vmath::Vadd(nPointsTot, Wzz, 1, tmp1, 1, tmp2, 1);
446  Vmath::Vmul(nPointsTot, m_zta[3], 1, tmp2, 1, tmp3, 1);
447  Vmath::Vadd(nPointsTot, Px, 1, tmp3, 1, Px, 1);
448 
449  Vmath::Smul(nPointsTot, m_kinvis, Px, 1, Px, 1); //* 1/Re
450 
451  //-------------------------------------------------------------------------
452  // y-component of the m_forcing - viscous component1:
453  // (V_z^'z^' - V_zz + ETA_tz^'z^')
454  //-------------------------------------------------------------------------
455  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
456  fields[1]->GetPhys(), tmp1); // Vz
457  fields[0]->HomogeneousBwdTrans(tmp1, tmp3); // Vz in physical space
458  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp3, tmp1); // Vzx
459  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], tmp3, tmp2); // Vzy
460  Vmath::Vmul(nPointsTot, m_zta[3], 1, tmp1, 1, tmp1, 1);
461  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp1, 1);
462  Vmath::Vmul(nPointsTot, m_eta[3], 1, tmp2, 1, tmp2, 1);
463  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp2, 1);
464  Vmath::Vsub(nPointsTot, Py, 1, tmp1, 1, Py, 1);
465  Vmath::Vsub(nPointsTot, Py, 1, tmp2, 1, Py, 1);
466 
467  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], V, tmp1); // Vx
468  Vmath::Vmul(nPointsTot, m_zta[4], 1, tmp1, 1, tmp2, 1);
469  Vmath::Vsub(nPointsTot, Py, 1, tmp2, 1, Py, 1);
470  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], tmp1, tmp2); // Vxx
471  Vmath::Vmul(nPointsTot, m_zta[9], 1, tmp2, 1, tmp3, 1);
472  Vmath::Vadd(nPointsTot, Py, 1, tmp3, 1, Py, 1);
473  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], tmp1, tmp2); // Vxy
474  Vmath::Vmul(nPointsTot, m_zta[8], 1, tmp2, 1, tmp3, 1);
475  Vmath::Smul(nPointsTot, 2.0, tmp3, 1, tmp3, 1);
476  Vmath::Vadd(nPointsTot, Py, 1, tmp3, 1, Py, 1);
477 
478  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], V, tmp1); // Vy
479  Vmath::Vmul(nPointsTot, m_eta[4], 1, tmp1, 1, tmp2, 1);
480  Vmath::Vsub(nPointsTot, Py, 1, tmp2, 1, Py, 1);
481  fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], tmp1, tmp2); // Vyy
482  Vmath::Vmul(nPointsTot, m_eta[9], 1, tmp2, 1, tmp3, 1);
483  Vmath::Vadd(nPointsTot, Py, 1, tmp3, 1, Py, 1);
484  Vmath::Vadd(nPointsTot, m_eta[7], 1, Py, 1, Py, 1);
485 
486  //-------------------------------------------------------------------------
487  // y-component of the m_forcing - viscous component2:
488  // ((ETA_z * W)_z^'z^' + ETA_z * (W_xx + W_yy))
489  //-------------------------------------------------------------------------
490  Vmath::Vmul(nPointsTot, m_eta[5], 1, W, 1, tmp1, 1);
491  Vmath::Vadd(nPointsTot, Py, 1, tmp1, 1, Py, 1);
492 
493  Vmath::Vmul(nPointsTot, m_zta[3], 1, Wx, 1, tmp1, 1);
494  Vmath::Vmul(nPointsTot, m_eta[4], 1, tmp1, 1, tmp2, 1);
495  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp3, 1);
496  Vmath::Vsub(nPointsTot, Py, 1, tmp3, 1, Py, 1);
497 
498  Vmath::Vmul(nPointsTot, m_zta[4], 1, Wx, 1, tmp1, 1);
499  Vmath::Vmul(nPointsTot, m_eta[3], 1, tmp1, 1, tmp2, 1);
500  Vmath::Vsub(nPointsTot, Py, 1, tmp2, 1, Py, 1);
501 
502  Vmath::Vmul(nPointsTot, m_eta[3], 1, Wy, 1, tmp1, 1);
503  Vmath::Vmul(nPointsTot, m_eta[4], 1, tmp1, 1, tmp2, 1);
504  Vmath::Smul(nPointsTot, 3.0, tmp2, 1, tmp3, 1);
505  Vmath::Vsub(nPointsTot, Py, 1, tmp3, 1, Py, 1);
506 
507  Vmath::Vmul(nPointsTot, m_eta[4], 1, Wz, 1, tmp1, 1);
508  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
509  Vmath::Vadd(nPointsTot, Py, 1, tmp2, 1, Py, 1);
510 
511  Vmath::Vmul(nPointsTot, m_eta[9], 1, Wxy, 1, tmp1, 1);
512  Vmath::Vmul(nPointsTot, m_zta[3], 1, tmp1, 1, tmp2, 1);
513  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp3, 1);
514  Vmath::Vadd(nPointsTot, Py, 1, tmp3, 1, Py, 1);
515 
516  Vmath::Vmul(nPointsTot, m_zta[8], 1, Wxz, 1, tmp1, 1);
517  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
518  Vmath::Vsub(nPointsTot, Py, 1, tmp2, 1, Py, 1);
519 
520  Vmath::Vmul(nPointsTot, m_eta[9], 1, Wyz, 1, tmp1, 1);
521  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp2, 1);
522  Vmath::Vsub(nPointsTot, Py, 1, tmp2, 1, Py, 1);
523 
524  Vmath::Vmul(nPointsTot, m_eta[3], 1, Wxx, 1, tmp1, 1);
525  Vmath::Vmul(nPointsTot, m_zta[9], 1, tmp1, 1, tmp2, 1);
526  Vmath::Vadd(nPointsTot, Py, 1, tmp2, 1, Py, 1);
527 
528  Vmath::Vmul(nPointsTot, m_eta[9], 1, m_eta[3], 1, tmp1, 1);
529  Vmath::Vmul(nPointsTot, tmp1, 1, Wyy, 1, tmp2, 1);
530  Vmath::Vadd(nPointsTot, Py, 1, tmp2, 1, Py, 1);
531 
532  Vmath::Vadd(nPointsTot, Wxx, 1, Wyy, 1, tmp1, 1);
533  Vmath::Vadd(nPointsTot, Wzz, 1, tmp1, 1, tmp2, 1);
534  Vmath::Vmul(nPointsTot, m_eta[3], 1, tmp2, 1, tmp3, 1);
535  Vmath::Vadd(nPointsTot, Py, 1, tmp3, 1, Py, 1);
536 
537  Vmath::Smul(nPointsTot, m_kinvis, Py, 1, Py, 1); //* 1/Re
538 
539  //-------------------------------------------------------------------------
540  // z-component of the m_forcing - viscous component: (W_z^'z^' - W_zz)
541  //-------------------------------------------------------------------------
542  Vmath::Vmul(nPointsTot, m_zta[3], 1, Wxz, 1, tmp1, 1);
543  Vmath::Smul(nPointsTot, 2.0, tmp1, 1, tmp1, 1);
544  Vmath::Vmul(nPointsTot, m_eta[3], 1, Wyz, 1, tmp2, 1);
545  Vmath::Smul(nPointsTot, 2.0, tmp2, 1, tmp2, 1);
546  Vmath::Vsub(nPointsTot, Pz, 1, tmp1, 1, Pz, 1);
547  Vmath::Vsub(nPointsTot, Pz, 1, tmp2, 1, Pz, 1);
548 
549  Vmath::Vmul(nPointsTot, m_zta[4], 1, Wx, 1, tmp2, 1);
550  Vmath::Vsub(nPointsTot, Pz, 1, tmp2, 1, Pz, 1);
551  Vmath::Vmul(nPointsTot, m_zta[9], 1, Wxx, 1, tmp3, 1);
552  Vmath::Vadd(nPointsTot, Pz, 1, tmp3, 1, Pz, 1);
553  Vmath::Vmul(nPointsTot, m_zta[8], 1, Wxy, 1, tmp3, 1);
554  Vmath::Smul(nPointsTot, 2.0, tmp3, 1, tmp3, 1);
555  Vmath::Vadd(nPointsTot, Pz, 1, tmp3, 1, Pz, 1);
556 
557  Vmath::Vmul(nPointsTot, m_eta[4], 1, Wy, 1, tmp2, 1);
558  Vmath::Vsub(nPointsTot, Pz, 1, tmp2, 1, Pz, 1);
559  Vmath::Vmul(nPointsTot, m_eta[9], 1, Wyy, 1, tmp3, 1);
560  Vmath::Vadd(nPointsTot, Pz, 1, tmp3, 1, Pz, 1);
561 
562  Vmath::Smul(nPointsTot, m_kinvis, Pz, 1, Pz, 1); //* 1/Re
563 
564  //-------------------------------------------------------------------------
565  // adding viscous and pressure components and transfroming back to wave
566  // space
567  //-------------------------------------------------------------------------
568  Vmath::Vadd(nPointsTot, Fx, 1, Px, 1, Fx, 1);
569  Vmath::Vadd(nPointsTot, Fy, 1, Py, 1, Fy, 1);
570  Vmath::Vadd(nPointsTot, Fz, 1, Pz, 1, Fz, 1);
571  fields[0]->HomogeneousFwdTrans(Fx, m_forcing[0]);
572  fields[0]->HomogeneousFwdTrans(Fy, m_forcing[1]);
573  fields[0]->HomogeneousFwdTrans(Fz, m_forcing[2]);
574 }
Nektar::ForcingMovingBody::m_kinvis
NekDouble m_kinvis
fluid viscous
Definition: ForcingMovingBody.h:138
Nektar::ForcingMovingBody::EvaluateAccelaration
void EvaluateAccelaration(const Array< OneD, NekDouble > &input, Array< OneD, NekDouble > &output, int npoints)
Definition: ForcingMovingBody.cpp:1593
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Nektar::ForcingMovingBody::m_eta
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
Definition: ForcingMovingBody.h:167
Nektar::ForcingMovingBody::m_forcing
Array< OneD, Array< OneD, NekDouble > > m_forcing
Definition: ForcingMovingBody.h:169
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::ForcingMovingBody::m_zta
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
Definition: ForcingMovingBody.h:165
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::ForcingMovingBody::CheckIsFromFile ( const TiXmlElement *  pForce)
private

Definition at line 1628 of file ForcingMovingBody.cpp.

References ASSERTL0, Nektar::LibUtilities::eFunctionTypeExpression, Nektar::LibUtilities::eFunctionTypeFile, m_funcName, m_IsFromFile, m_motion, and Nektar::SolverUtils::Forcing::m_session.

Referenced by v_InitObject().

1629 {
1630 
1631  m_funcName = Array<OneD, std::string> (3);
1632  m_motion = Array<OneD, std::string> (2);
1633  m_motion[0] = "x";
1634  m_motion[1] = "y";
1635 
1636  m_IsFromFile = Array<OneD, bool> (6);
1637  // Loading the x-dispalcement (m_zta) and the y-displacement (m_eta)
1638  // Those two variables are bith functions of z and t and the may come
1639  // from an equation (forced vibration) or from another solver which, given
1640  // the aerodynamic forces at the previous step, calculates the
1641  // displacements.
1642 
1643  //Get the body displacement: m_zta and m_eta
1644  const TiXmlElement* funcNameElmt_D
1645  = pForce->FirstChildElement("DISPLACEMENTS");
1646  ASSERTL0(funcNameElmt_D,
1647  "MOVINGBODYFORCE tag has to specify a function name which "
1648  "prescribes the body displacement as d(z,t).");
1649 
1650  m_funcName[0] = funcNameElmt_D->GetText();
1651  ASSERTL0(m_session->DefinesFunction(m_funcName[0]),
1652  "Function '" + m_funcName[0] + "' not defined.");
1653 
1654  //Get the body velocity of movement: d(m_zta)/dt and d(m_eta)/dt
1655  const TiXmlElement* funcNameElmt_V
1656  = pForce->FirstChildElement("VELOCITIES");
1657  ASSERTL0(funcNameElmt_D,
1658  "MOVINGBODYFORCE tag has to specify a function name which "
1659  "prescribes the body velocity of movement as v(z,t).");
1660 
1661  m_funcName[1] = funcNameElmt_V->GetText();
1662  ASSERTL0(m_session->DefinesFunction(m_funcName[1]),
1663  "Function '" + m_funcName[1] + "' not defined.");
1664 
1665 
1666  //Get the body acceleration: dd(m_zta)/ddt and dd(m_eta)/ddt
1667  const TiXmlElement* funcNameElmt_A
1668  = pForce->FirstChildElement("ACCELERATIONS");
1669  ASSERTL0(funcNameElmt_A,
1670  "MOVINGBODYFORCE tag has to specify a function name which "
1671  "prescribes the body acceleration as a(z,t).");
1672 
1673  m_funcName[2] = funcNameElmt_A->GetText();
1674  ASSERTL0(m_session->DefinesFunction(m_funcName[2]),
1675  "Function '" + m_funcName[2] + "' not defined.");
1676 
1677  LibUtilities::FunctionType vType;
1678 
1679  // Check Displacement x
1680  vType = m_session->GetFunctionType(m_funcName[0],m_motion[0]);
1681  if(vType == LibUtilities::eFunctionTypeExpression)
1682  {
1683  m_IsFromFile[0] = false;
1684  }
1685  else if (vType == LibUtilities::eFunctionTypeFile)
1686  {
1687  m_IsFromFile[0] = true;
1688  }
1689  else
1690  {
1691  ASSERTL0(false, "The displacements in x must be specified via an "
1692  "equation <E> or a file <F>");
1693  }
1694 
1695  // Check Displacement y
1696  vType = m_session->GetFunctionType(m_funcName[0],m_motion[1]);
1697  if(vType == LibUtilities::eFunctionTypeExpression)
1698  {
1699  m_IsFromFile[1] = false;
1700  }
1701  else if (vType == LibUtilities::eFunctionTypeFile)
1702  {
1703  m_IsFromFile[1] = true;
1704  }
1705  else
1706  {
1707  ASSERTL0(false, "The displacements in y must be specified via an "
1708  "equation <E> or a file <F>");
1709  }
1710 
1711  // Check Velocity x
1712  vType = m_session->GetFunctionType(m_funcName[1],m_motion[0]);
1713  if (vType == LibUtilities::eFunctionTypeExpression)
1714  {
1715  m_IsFromFile[2] = false;
1716  }
1717  else if (vType == LibUtilities::eFunctionTypeFile)
1718  {
1719  m_IsFromFile[2] = true;
1720  }
1721  else
1722  {
1723  ASSERTL0(false, "The velocities in x must be specified via an equation "
1724  "<E> or a file <F>");
1725  }
1726 
1727  // Check Velocity y
1728  vType = m_session->GetFunctionType(m_funcName[1],m_motion[1]);
1729  if (vType == LibUtilities::eFunctionTypeExpression)
1730  {
1731  m_IsFromFile[3] = false;
1732  }
1733  else if (vType == LibUtilities::eFunctionTypeFile)
1734  {
1735  m_IsFromFile[3] = true;
1736  }
1737  else
1738  {
1739  ASSERTL0(false, "The velocities in y must be specified via an equation "
1740  "<E> or a file <F>");
1741  }
1742 
1743  // Check Acceleration x
1744  vType = m_session->GetFunctionType(m_funcName[2],m_motion[0]);
1745  if (vType == LibUtilities::eFunctionTypeExpression)
1746  {
1747  m_IsFromFile[4] = false;
1748  }
1749  else if (vType == LibUtilities::eFunctionTypeFile)
1750  {
1751  m_IsFromFile[4] = true;
1752  }
1753  else
1754  {
1755  ASSERTL0(false, "The accelerations in x must be specified via an "
1756  "equation <E> or a file <F>");
1757  }
1758 
1759  // Check Acceleration y
1760  vType = m_session->GetFunctionType(m_funcName[2],m_motion[1]);
1761  if (vType == LibUtilities::eFunctionTypeExpression)
1762  {
1763  m_IsFromFile[5] = false;
1764  }
1765  else if (vType == LibUtilities::eFunctionTypeFile)
1766  {
1767  m_IsFromFile[5] = true;
1768  }
1769  else
1770  {
1771  ASSERTL0(false, "The accelerations in y must be specified via an "
1772  "equation <E> or a file <F>");
1773  }
1774 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ForcingMovingBody::m_motion
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
Definition: ForcingMovingBody.h:161
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::ForcingMovingBody::m_IsFromFile
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
Definition: ForcingMovingBody.h:163
Nektar::ForcingMovingBody::m_funcName
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
Definition: ForcingMovingBody.h:159
static SolverUtils::ForcingSharedPtr Nektar::ForcingMovingBody::create ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const unsigned int &  pNumForcingFields,
const TiXmlElement *  pForce 
)
inlinestatic

Creates an instance of this class.

Definition at line 56 of file ForcingMovingBody.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and Nektar::SolverUtils::ForcingSharedPtr.

61  {
62  SolverUtils::ForcingSharedPtr p =
63  MemoryManager<ForcingMovingBody>::
64  AllocateSharedPtr(pSession);
65  p->InitObject(pFields, pNumForcingFields, pForce);
66  return p;
67  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::SolverUtils::ForcingSharedPtr
SOLVER_UTILS_EXPORT typedef boost::shared_ptr< Forcing > ForcingSharedPtr
A shared pointer to an EquationSystem object.
Definition: Forcing.h:51
void Nektar::ForcingMovingBody::EvaluateAccelaration ( const Array< OneD, NekDouble > &  input,
Array< OneD, NekDouble > &  output,
int  npoints 
)
private

Definition at line 1593 of file ForcingMovingBody.cpp.

References m_movingBodyCalls, m_timestep, m_W, RollOver(), Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

Referenced by CalculateForcing().

1597 {
1598  m_movingBodyCalls++;
1599 
1600  // Rotate acceleration term
1601  RollOver(m_W);
1602 
1603  Vmath::Vcopy(npoints, input, 1, m_W[0], 1);
1604 
1605  //Calculate acceleration term at level n based on previous steps
1606  if (m_movingBodyCalls > 1)
1607  {
1608  Vmath::Smul(npoints,
1609  1.0,
1610  m_W[0], 1,
1611  output, 1);
1612  Vmath::Svtvp(npoints,
1613  -1.0,
1614  m_W[1], 1,
1615  output, 1,
1616  output, 1);
1617  Vmath::Smul(npoints,
1618  1.0/m_timestep,
1619  output, 1,
1620  output, 1);
1621  }
1622 }
Nektar::ForcingMovingBody::m_W
Array< OneD, Array< OneD, NekDouble > > m_W
srorage for the velocity in z-direction
Definition: ForcingMovingBody.h:149
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Nektar::ForcingMovingBody::m_movingBodyCalls
int m_movingBodyCalls
number of times the movbody have been called
Definition: ForcingMovingBody.h:131
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::ForcingMovingBody::m_timestep
NekDouble m_timestep
time step
Definition: ForcingMovingBody.h:139
Nektar::ForcingMovingBody::RollOver
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
Definition: ForcingMovingBody.cpp:1574
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void Nektar::ForcingMovingBody::EvaluateStructDynModel ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
NekDouble  time 
)
private

Definition at line 921 of file ForcingMovingBody.cpp.

References Vmath::Fill(), m_Aeroforces, m_eta, m_fA, m_fV, m_motion, m_MotionVars, m_movingBodyCalls, m_np, Nektar::SolverUtils::Forcing::m_session, m_vdim, m_zta, RollOver(), Vmath::Smul(), Vmath::Svtvp(), TensionedCableModel(), and Vmath::Vcopy().

Referenced by UpdateMotion().

924 {
925  //Get the hydrodynamic forces from the fluid solver
926  Array<OneD, Array <OneD, NekDouble> > fces(m_motion.num_elements());
927 
928  fces[0] = Array <OneD, NekDouble> (m_np,0.0);
929  fces[1] = Array <OneD, NekDouble> (m_np,0.0);
930 
931  for(int plane = 0; plane < m_np; plane++)
932  {
933  fces[0][plane] = m_Aeroforces[plane];
934  fces[1][plane] = m_Aeroforces[plane+m_np];
935  }
936 
937  // Fictitious mass method used to stablize the explicit coupling at
938  // relatively lower mass ratio
939  bool fictmass;
940  m_session->MatchSolverInfo("FictitiousMassMethod", "True",
941  fictmass, false);
942  if(fictmass)
943  {
944  NekDouble fictrho, fictdamp;
945  m_session->LoadParameter("FictMass", fictrho);
946  m_session->LoadParameter("FictDamp", fictdamp);
947 
948  static NekDouble Betaq_Coeffs[2][2] =
949  {{1.0, 0.0},{2.0, -1.0}};
950 
951  // only consider second order approximation for fictitious variables
952  int intOrder= 2;
953  int nint = min(m_movingBodyCalls+1,intOrder);
954  int nlevels = m_fV[0].num_elements();
955 
956  for(int i = 0; i < m_motion.num_elements(); ++i)
957  {
958  RollOver(m_fV[i]);
959  RollOver(m_fA[i]);
960 
961  int Voffset = i*3*m_np+m_np;
962  int Aoffset = i*3*m_np+2*m_np;
963 
964  Vmath::Vcopy(m_np, m_MotionVars+Voffset, 1, m_fV[i][0], 1);
965  Vmath::Vcopy(m_np, m_MotionVars+Aoffset, 1, m_fA[i][0], 1);
966 
967  // Extrapolate to n+1
968  Vmath::Smul(m_np,
969  Betaq_Coeffs[nint-1][nint-1],
970  m_fV[i][nint-1], 1,
971  m_fV[i][nlevels-1], 1);
972  Vmath::Smul(m_np,
973  Betaq_Coeffs[nint-1][nint-1],
974  m_fA[i][nint-1], 1,
975  m_fA[i][nlevels-1], 1);
976 
977  for(int n = 0; n < nint-1; ++n)
978  {
979  Vmath::Svtvp(m_np,
980  Betaq_Coeffs[nint-1][n],
981  m_fV[i][n],1,m_fV[i][nlevels-1],1,
982  m_fV[i][nlevels-1],1);
983  Vmath::Svtvp(m_np,
984  Betaq_Coeffs[nint-1][n],
985  m_fA[i][n],1,m_fA[i][nlevels-1],1,
986  m_fA[i][nlevels-1],1);
987  }
988 
989  // Add the fictitious forces on the RHS of the equation
990  Vmath::Svtvp(m_np, fictdamp,m_fV[i][nlevels-1],1,
991  fces[i],1,fces[i],1);
992  Vmath::Svtvp(m_np, fictrho, m_fA[i][nlevels-1],1,
993  fces[i],1,fces[i],1);
994  }
995  }
996 
997  // Tensioned cable model is evaluated in wave space
998  for(int n = 0, cn = 1; n < m_vdim; n++, cn--)
999  {
1000  int offset = cn*3*m_np;
1001  Array<OneD, NekDouble> tmp(3*m_np);
1002  TensionedCableModel(pFields, fces[cn],
1003  tmp = m_MotionVars+offset);
1004  }
1005 
1006  // Set the m_forcing term based on the motion of the cable
1007  for(int var = 0; var < 3; var++)
1008  {
1009  for(int plane = 0; plane < m_np; plane++)
1010  {
1011  int n = pFields[0]->GetPlane(plane)->GetTotPoints();
1012 
1013  Array<OneD, NekDouble> tmp;
1014 
1015  int offset = plane * n;
1016  int xoffset = var * m_np+plane;
1017  int yoffset = 3*m_np + xoffset;
1018 
1019  Vmath::Fill(n, m_MotionVars[xoffset], tmp = m_zta[var] + offset, 1);
1020  Vmath::Fill(n, m_MotionVars[yoffset], tmp = m_eta[var] + offset, 1);
1021  }
1022  }
1023 }
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Nektar::ForcingMovingBody::m_Aeroforces
Array< OneD, NekDouble > m_Aeroforces
storage for the cable's force(x,y) variables
Definition: ForcingMovingBody.h:145
Nektar::ForcingMovingBody::m_movingBodyCalls
int m_movingBodyCalls
number of times the movbody have been called
Definition: ForcingMovingBody.h:131
Nektar::ForcingMovingBody::m_motion
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
Definition: ForcingMovingBody.h:161
Nektar::ForcingMovingBody::m_MotionVars
Array< OneD, NekDouble > m_MotionVars
storage for the cable's motion(x,y) variables
Definition: ForcingMovingBody.h:147
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::ForcingMovingBody::m_vdim
int m_vdim
vibration dimension
Definition: ForcingMovingBody.h:133
Nektar::ForcingMovingBody::TensionedCableModel
void TensionedCableModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
Definition: ForcingMovingBody.cpp:580
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::ForcingMovingBody::m_eta
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
Definition: ForcingMovingBody.h:167
Nektar::ForcingMovingBody::m_fV
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
Definition: ForcingMovingBody.h:151
Nektar::ForcingMovingBody::m_np
int m_np
number of planes per processors
Definition: ForcingMovingBody.h:132
Nektar::ForcingMovingBody::m_fA
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
Definition: ForcingMovingBody.h:153
Nektar::ForcingMovingBody::RollOver
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
Definition: ForcingMovingBody.cpp:1574
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Nektar::ForcingMovingBody::m_zta
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
Definition: ForcingMovingBody.h:165
void Nektar::ForcingMovingBody::InitialiseCableModel ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields 
)
private

Definition at line 1106 of file ForcingMovingBody.cpp.

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::LibUtilities::GetNektarFFTFactory(), m_Aeroforces, m_fA, m_FFT, m_funcName, m_fV, m_IsFromFile, m_kinvis, m_lhom, m_motion, m_MotionVars, m_movingBodyCalls, m_np, Nektar::SolverUtils::Forcing::m_session, m_structdamp, m_structrho, m_timestep, m_vdim, m_W, Vmath::Sadd(), SetDynEqCoeffMatrix(), and Vmath::Smul().

Referenced by v_InitObject().

1109 {
1110  // storage of spanwise-velocity for current and previous time levels
1111  // used to calculate dw/dt in m_forcing term
1112  int nPointsTot = pFields[0]->GetNpoints();
1113  m_W = Array<OneD, Array<OneD, NekDouble> > (2);
1114  for(int n = 0; n < 2; ++n)
1115  {
1116  m_W[n] = Array<OneD, NekDouble> (nPointsTot, 0.0);
1117  }
1118 
1119  m_session->LoadParameter("Kinvis",m_kinvis);
1120  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
1121 
1122  m_movingBodyCalls = 0;
1123 
1124  Array<OneD, unsigned int> ZIDs;
1125  ZIDs = pFields[0]->GetZIDs();
1126  m_np = ZIDs.num_elements();
1127 
1128  m_Aeroforces = Array<OneD, NekDouble>(2*m_np,0.0);
1129  m_MotionVars = Array<OneD, NekDouble>(6*m_np,0.0);
1130 
1131  std::string vibtype = m_session->GetSolverInfo("VibrationType");
1132 
1133  if(boost::iequals(vibtype, "Constrained"))
1134  {
1135  m_vdim = 1;
1136  }
1137  else if (boost::iequals(vibtype, "Free"))
1138  {
1139  m_vdim = 2;
1140  }
1141  else if (boost::iequals(vibtype, "Forced"))
1142  {
1143  return;
1144  }
1145 
1146  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
1147 
1148  bool homostrip;
1149  m_session->MatchSolverInfo("HomoStrip","True",homostrip,false);
1150 
1151  if(!homostrip)
1152  {
1153  m_session->LoadParameter("LZ", m_lhom);
1154  int nplanes = m_session->GetParameter("HomModesZ");
1155  m_FFT =
1156  LibUtilities::GetNektarFFTFactory().CreateInstance(
1157  "NekFFTW", nplanes);
1158  }
1159  else
1160  {
1161  int nstrips;
1162 
1163  NekDouble DistStrip;
1164 
1165  m_session->LoadParameter("DistStrip", DistStrip);
1166  m_session->LoadParameter("Strip_Z", nstrips);
1167  m_lhom = nstrips * DistStrip;
1168  m_FFT =
1169  LibUtilities::GetNektarFFTFactory().CreateInstance(
1170  "NekFFTW", nstrips);
1171  }
1172 
1173  // load the structural dynamic parameters from xml file
1174  m_session->LoadParameter("StructRho", m_structrho);
1175  m_session->LoadParameter("StructDamp", m_structdamp, 0.0);
1176 
1177  // Identify whether the fictitious mass method is active for explicit
1178  // coupling of fluid solver and structural dynamics solver
1179  bool fictmass;
1180  m_session->MatchSolverInfo("FictitiousMassMethod", "True",
1181  fictmass, false);
1182  if(fictmass)
1183  {
1184  NekDouble fictrho, fictdamp;
1185  m_session->LoadParameter("FictMass", fictrho);
1186  m_session->LoadParameter("FictDamp", fictdamp);
1187  m_structrho += fictrho;
1188  m_structdamp += fictdamp;
1189 
1190  // Storage array of Struct Velocity and Acceleration used for
1191  // extrapolation of fictitious force
1192  m_fV = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
1193  m_motion.num_elements());
1194  m_fA = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
1195  m_motion.num_elements());
1196  for (int i = 0; i < m_motion.num_elements(); ++i)
1197  {
1198  m_fV[i] = Array<OneD, Array<OneD, NekDouble> > (2);
1199  m_fA[i] = Array<OneD, Array<OneD, NekDouble> > (2);
1200 
1201  for(int n = 0; n < 2; ++n)
1202  {
1203  m_fV[i][n] = Array<OneD, NekDouble>(m_np, 0.0);
1204  m_fA[i][n] = Array<OneD, NekDouble>(m_np, 0.0);
1205  }
1206  }
1207  }
1208 
1209  // Setting the coefficient matrices for solving structural dynamic ODEs
1210  SetDynEqCoeffMatrix(pFields);
1211 
1212  // Set initial condition for cable's motion
1213  int cnt = 0;
1214 
1215  for(int j = 0; j < m_funcName.num_elements(); j++)
1216  {
1217 
1218  // loading from the specified files through inputstream
1219  if(m_IsFromFile[cnt] && m_IsFromFile[cnt+1])
1220  {
1221 
1222  std::ifstream inputStream;
1223 
1224  int nzpoints;
1225 
1226  if(homostrip)
1227  {
1228  int nstrips;
1229  m_session->LoadParameter("Strip_Z", nstrips);
1230  nzpoints = nstrips;
1231  }
1232  else
1233  {
1234  nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
1235  }
1236 
1237 
1238  if (vcomm->GetRank() == 0)
1239  {
1240 
1241  Array<OneD, NekDouble> motion_x(3*nzpoints, 0.0);
1242  Array<OneD, NekDouble> motion_y(3*nzpoints, 0.0);
1243 
1244  Array<OneD, std::string> tmp(9);
1245 
1246  std::string filename = m_session->GetFunctionFilename(
1247  m_funcName[j], m_motion[0]);
1248 
1249  // Open intputstream for cable motions
1250  inputStream.open(filename.c_str());
1251 
1252  // Import the head string from the file
1253  for(int n = 0; n < tmp.num_elements(); n++)
1254  {
1255  inputStream >> tmp[n];
1256  }
1257 
1258  NekDouble time, z_cds;
1259  // Import the motion variables from the file
1260  for (int n = 0; n < nzpoints; n++)
1261  {
1262  inputStream >> setprecision(6) >> time;
1263  inputStream >> setprecision(6) >> z_cds;
1264  inputStream >> setprecision(8) >> motion_x[n];
1265  inputStream >> setprecision(8) >> motion_x[n+nzpoints];
1266  inputStream >> setprecision(8) >> motion_x[n+2*nzpoints];
1267  inputStream >> setprecision(8) >> motion_y[n];
1268  inputStream >> setprecision(8) >> motion_y[n+nzpoints];
1269  inputStream >> setprecision(8) >> motion_y[n+2*nzpoints];
1270  }
1271 
1272  // Close inputstream for cable motions
1273  inputStream.close();
1274 
1275  if(!homostrip)
1276  {
1277  int nproc = vcomm->GetColumnComm()->GetSize();
1278  for (int i = 1; i < nproc; ++i)
1279  {
1280  for(int plane = 0; plane < m_np; plane++)
1281  {
1282  for (int var = 0; var < 3; var++)
1283  {
1284  int xoffset = var*m_np+plane;
1285  int yoffset = 3*m_np+xoffset;
1286  m_MotionVars[xoffset] =
1287  motion_x[plane+i*m_np+var*nzpoints];
1288  m_MotionVars[yoffset] =
1289  motion_y[plane+i*m_np+var*nzpoints];
1290  }
1291  }
1292 
1293  vcomm->GetColumnComm()->Send(i, m_MotionVars);
1294  }
1295 
1296  for(int plane = 0; plane < m_np; plane++)
1297  {
1298  for (int var = 0; var < 3; var++)
1299  {
1300  int xoffset = var*m_np+plane;
1301  int yoffset = 3*m_np+xoffset;
1302  m_MotionVars[xoffset] = motion_x[plane+var*nzpoints];
1303  m_MotionVars[yoffset] = motion_y[plane+var*nzpoints];
1304  }
1305  }
1306  }
1307  else //for strip model, make sure that each processor gets its motion values correctly
1308  {
1309  int nstrips;
1310  m_session->LoadParameter("Strip_Z", nstrips);
1311 
1312  int nproc = vcomm->GetColumnComm()->GetSize();
1313 
1314  for (int i = 1; i < nstrips; ++i)
1315  {
1316  for(int plane = 0; plane < m_np; plane++)
1317  {
1318  for (int var = 0; var < 3; var++)
1319  {
1320  int xoffset = var*m_np+plane;
1321  int yoffset = 3*m_np+xoffset;
1322  m_MotionVars[xoffset] = motion_x[i+var*nstrips];
1323  m_MotionVars[yoffset] = motion_y[i+var*nstrips];
1324  }
1325  }
1326 
1327  for (int j = 0; j < nproc/nstrips; j++)
1328  {
1329  vcomm->GetColumnComm()->Send(i+j*nstrips, m_MotionVars);
1330  }
1331  }
1332 
1333  for(int plane = 0; plane < m_np; plane++)
1334  {
1335  for (int var = 0; var < 3; var++)
1336  {
1337  int xoffset = var*m_np+plane;
1338  int yoffset = 3*m_np+xoffset;
1339  m_MotionVars[xoffset] = motion_x[var*nstrips];
1340  m_MotionVars[yoffset] = motion_y[var*nstrips];
1341  }
1342  }
1343 
1344  for (int j = 1; j < nproc/nstrips; j++)
1345  {
1346  vcomm->GetColumnComm()->Send(j*nstrips, m_MotionVars);
1347  }
1348  }
1349  }
1350  else
1351  {
1352  if(!homostrip)
1353  {
1354  int colrank = vcomm->GetColumnComm()->GetRank();
1355  int nproc = vcomm->GetColumnComm()->GetSize();
1356 
1357  for (int j = 1; j < nproc; j++)
1358  {
1359  if(colrank == j)
1360  {
1361  vcomm->GetColumnComm()->Recv(0, m_MotionVars);
1362  }
1363  }
1364  }
1365  else //for strip model
1366  {
1367  int nstrips;
1368  m_session->LoadParameter("Strip_Z", nstrips);
1369 
1370  int colrank = vcomm->GetColumnComm()->GetRank();
1371  int nproc = vcomm->GetColumnComm()->GetSize();
1372 
1373  for(int i = 1; i < nstrips; i++)
1374  {
1375  for (int j = 0; j < nproc/nstrips; j++)
1376  {
1377  if(colrank == i+j*nstrips)
1378  {
1379  vcomm->GetColumnComm()->Recv(0, m_MotionVars);
1380  }
1381  }
1382  }
1383 
1384  for (int j = 1; j < nproc/nstrips; j++)
1385  {
1386  if(colrank == j*nstrips)
1387  {
1388  vcomm->GetColumnComm()->Recv(0, m_MotionVars);
1389  }
1390  }
1391  }
1392  }
1393 
1394  cnt = cnt + 2;
1395  }
1396  else //Evaluate from the functions specified in xml file
1397  {
1398  Array<OneD, NekDouble> x0(m_np, 0.0);
1399  Array<OneD, NekDouble> x1(m_np, 0.0);
1400  Array<OneD, NekDouble> x2(m_np, 0.0);
1401 
1402  if(!homostrip)
1403  {
1404  int nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
1405  Array<OneD, NekDouble> z_coords(nzpoints,0.0);
1406  Array<OneD, const NekDouble> pts
1407  = pFields[0]->GetHomogeneousBasis()->GetZ();
1408 
1409  Vmath::Smul(nzpoints,m_lhom/2.0,pts,1,z_coords,1);
1410  Vmath::Sadd(nzpoints,m_lhom/2.0,z_coords,1,z_coords,1);
1411 
1412  for (int plane = 0; plane < m_np; plane++)
1413  {
1414  x2[plane] = z_coords[ZIDs[plane]];
1415  }
1416  }
1417  else
1418  {
1419  int nstrips;
1420  m_session->LoadParameter("Strip_Z", nstrips);
1421 
1422  ASSERTL0(m_session->DefinesSolverInfo("USEFFT"),
1423  "Fourier transformation of cable motion is currently "
1424  "implemented only for FFTW module.");
1425 
1426  NekDouble DistStrip;
1427  m_session->LoadParameter("DistStrip", DistStrip);
1428 
1429  int colrank = vcomm->GetColumnComm()->GetRank();
1430  int nproc = vcomm->GetColumnComm()->GetSize();
1431 
1432  for(int i = 0; i < nstrips; ++i)
1433  {
1434  for(int j = 0; j < nproc/nstrips; j++)
1435  {
1436  if (colrank == i+j*nstrips)
1437  {
1438  for (int plane = 0; plane < m_np; plane++)
1439  {
1440  x2[plane] = i*DistStrip;
1441  }
1442  }
1443  }
1444  }
1445  }
1446 
1447  int xoffset = j*m_np;
1448  int yoffset = 3*m_np+xoffset;
1449 
1450  Array<OneD, NekDouble> tmp(m_np,0.0);
1451  LibUtilities::EquationSharedPtr ffunc0,ffunc1;
1452 
1453  ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
1454  ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
1455 
1456  ffunc0->Evaluate(x0, x1, x2, 0.0, tmp = m_MotionVars+xoffset);
1457  ffunc1->Evaluate(x0, x1, x2, 0.0, tmp = m_MotionVars+yoffset);
1458  cnt = cnt + 2;
1459  }
1460  }
1461 }
Nektar::ForcingMovingBody::m_structrho
NekDouble m_structrho
mass of the cable per unit length
Definition: ForcingMovingBody.h:135
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ForcingMovingBody::m_structdamp
NekDouble m_structdamp
damping ratio of the cable
Definition: ForcingMovingBody.h:136
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Nektar::ForcingMovingBody::SetDynEqCoeffMatrix
void SetDynEqCoeffMatrix(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Definition: ForcingMovingBody.cpp:1467
Nektar::ForcingMovingBody::m_W
Array< OneD, Array< OneD, NekDouble > > m_W
srorage for the velocity in z-direction
Definition: ForcingMovingBody.h:149
Nektar::ForcingMovingBody::m_Aeroforces
Array< OneD, NekDouble > m_Aeroforces
storage for the cable's force(x,y) variables
Definition: ForcingMovingBody.h:145
Nektar::LibUtilities::GetNektarFFTFactory
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
Nektar::ForcingMovingBody::m_movingBodyCalls
int m_movingBodyCalls
number of times the movbody have been called
Definition: ForcingMovingBody.h:131
Nektar::ForcingMovingBody::m_kinvis
NekDouble m_kinvis
fluid viscous
Definition: ForcingMovingBody.h:138
Nektar::ForcingMovingBody::m_motion
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
Definition: ForcingMovingBody.h:161
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::ForcingMovingBody::m_MotionVars
Array< OneD, NekDouble > m_MotionVars
storage for the cable's motion(x,y) variables
Definition: ForcingMovingBody.h:147
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::ForcingMovingBody::m_vdim
int m_vdim
vibration dimension
Definition: ForcingMovingBody.h:133
Nektar::ForcingMovingBody::m_FFT
LibUtilities::NektarFFTSharedPtr m_FFT
Definition: ForcingMovingBody.h:141
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::ForcingMovingBody::m_timestep
NekDouble m_timestep
time step
Definition: ForcingMovingBody.h:139
Nektar::LibUtilities::EquationSharedPtr
boost::shared_ptr< Equation > EquationSharedPtr
Definition: SessionReader.h:100
Nektar::ForcingMovingBody::m_IsFromFile
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
Definition: ForcingMovingBody.h:163
Nektar::ForcingMovingBody::m_lhom
NekDouble m_lhom
length ratio of the cable
Definition: ForcingMovingBody.h:137
Nektar::ForcingMovingBody::m_fV
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
Definition: ForcingMovingBody.h:151
Nektar::ForcingMovingBody::m_np
int m_np
number of planes per processors
Definition: ForcingMovingBody.h:132
Nektar::ForcingMovingBody::m_funcName
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
Definition: ForcingMovingBody.h:159
Nektar::ForcingMovingBody::m_fA
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
Definition: ForcingMovingBody.h:153
void Nektar::ForcingMovingBody::InitialiseFilter ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const TiXmlElement *  pForce 
)
private

Definition at line 1779 of file ForcingMovingBody.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, and m_MovBodyfilter.

Referenced by v_InitObject().

1783 {
1784  // Get the outputfile name, output frequency and
1785  // the boundary's ID for the cable's wall
1786  std::string typeStr = pForce->Attribute("TYPE");
1787  std::map<std::string, std::string> vParams;
1788 
1789  const TiXmlElement *param = pForce->FirstChildElement("PARAM");
1790  while (param)
1791  {
1792  ASSERTL0(param->Attribute("NAME"),
1793  "Missing attribute 'NAME' for parameter in filter "
1794  + typeStr + "'.");
1795  std::string nameStr = param->Attribute("NAME");
1796 
1797  ASSERTL0(param->GetText(), "Empty value string for param.");
1798  std::string valueStr = param->GetText();
1799 
1800  vParams[nameStr] = valueStr;
1801 
1802  param = param->NextSiblingElement("PARAM");
1803  }
1804 
1805  // Creat the filter for MovingBody, where we performed the calculation of
1806  // fluid forces and write both the aerodynamic forces and motion variables
1807  // into the output files
1808  m_MovBodyfilter = MemoryManager<FilterMovingBody>::
1809  AllocateSharedPtr(pSession, vParams);
1810 
1811  // Initialise the object of MovingBody filter
1812  m_MovBodyfilter->Initialise(pFields, 0.0);
1813 
1814 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::ForcingMovingBody::m_MovBodyfilter
FilterMovingBodySharedPtr m_MovBodyfilter
Definition: ForcingMovingBody.h:143
void Nektar::ForcingMovingBody::MappingBndConditions ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pfields,
const Array< OneD, Array< OneD, NekDouble > > &  fields,
NekDouble  time 
)
private

Definition at line 1028 of file ForcingMovingBody.cpp.

References Nektar::LibUtilities::Equation::Evaluate(), Vmath::Fill(), Nektar::NekConstants::kNekUnsetDouble, m_lhom, m_motion, m_MotionVars, m_np, Vmath::Vcopy(), and Vmath::Vsub().

Referenced by v_Apply().

1032 {
1033  Array<OneD, Array<OneD, NekDouble> > BndV (
1034  m_motion.num_elements());
1035  for (int i = 0; i < m_motion.num_elements(); ++i)
1036  {
1037  BndV[i] = Array<OneD, NekDouble>(m_np, 0.0);
1038  }
1039 
1040  for(int i = 0; i < m_motion.num_elements(); ++i)
1041  {
1042  int offset = i*3*m_np+m_np;
1043  Vmath::Vcopy(m_np, m_MotionVars+offset, 1, BndV[i], 1);
1044  }
1045 
1046  Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExps;
1047  Array<OneD, const SpatialDomains::BoundaryConditionShPtr> bndConds;
1048 
1049  const Array<OneD, const NekDouble> z
1050  = pFields[0]->GetHomogeneousBasis()->GetZ();
1051  int nbnd = pFields[0]->GetBndCondExpansions().num_elements();
1052 
1053  for (int i = 0; i < nbnd; ++i)
1054  {
1055  for (int plane = 0; plane < m_np; plane++)
1056  {
1057  for ( int dim = 0; dim < m_motion.num_elements(); dim++)
1058  {
1059  bndCondExps = pFields[dim]->GetPlane(plane)
1060  ->GetBndCondExpansions();
1061  bndConds = pFields[dim]->GetPlane(plane)->GetBndConditions();
1062  if (boost::iequals(bndConds[i]->GetUserDefined(),"MovingBody"))
1063  {
1064  int npoints = bndCondExps[i]->GetNpoints();
1065  Array<OneD, NekDouble> x0(npoints, 0.0);
1066  Array<OneD, NekDouble> x1(npoints, 0.0);
1067  Array<OneD, NekDouble> x2(npoints, 0.0);
1068  Array<OneD, NekDouble> tmp(npoints,0.0);
1069 
1070  NekDouble local_z = z[pFields[0]->GetTransposition()
1071  ->GetPlaneID(plane)];
1072  NekDouble x2_in = 0.5*m_lhom*(1.0+local_z);
1073  // Homogeneous input case for x2.
1074  if (x2_in == NekConstants::kNekUnsetDouble)
1075  {
1076  bndCondExps[i]->GetCoords(x0,x1,x2);
1077  }
1078  else
1079  {
1080  bndCondExps[i]->GetCoords(x0, x1, x2);
1081  Vmath::Fill(npoints, x2_in, x2, 1);
1082  }
1083 
1084  LibUtilities::Equation condition =
1085  boost::static_pointer_cast<
1086  SpatialDomains::DirichletBoundaryCondition>
1087  (bndConds[i])->
1088  m_dirichletCondition;
1089 
1090  condition.Evaluate(x0, x1, x2, time,
1091  bndCondExps[i]->UpdatePhys());
1092  Vmath::Fill(npoints, BndV[dim][plane], tmp, 1);
1093  Vmath::Vsub(npoints, bndCondExps[i]->UpdatePhys(), 1,
1094  tmp, 1,
1095  bndCondExps[i]->UpdatePhys(), 1);
1096  }
1097  }
1098  }
1099  }
1100 }
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Nektar::ForcingMovingBody::m_motion
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
Definition: ForcingMovingBody.h:161
Nektar::ForcingMovingBody::m_MotionVars
Array< OneD, NekDouble > m_MotionVars
storage for the cable's motion(x,y) variables
Definition: ForcingMovingBody.h:147
Nektar::LibUtilities::Equation::Evaluate
NekDouble Evaluate() const
Definition: Equation.h:102
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition: NektarUnivConsts.hpp:47
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Nektar::ForcingMovingBody::m_lhom
NekDouble m_lhom
length ratio of the cable
Definition: ForcingMovingBody.h:137
Nektar::ForcingMovingBody::m_np
int m_np
number of planes per processors
Definition: ForcingMovingBody.h:132
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void Nektar::ForcingMovingBody::RollOver ( Array< OneD, Array< OneD, NekDouble > > &  input)
private

Function to roll time-level storages to the next step layout. The stored data associated with the oldest time-level (not required anymore) are moved to the top, where they will be overwritten as the solution process progresses.

Definition at line 1574 of file ForcingMovingBody.cpp.

Referenced by EvaluateAccelaration(), and EvaluateStructDynModel().

1575 {
1576  int nlevels = input.num_elements();
1577 
1578  Array<OneD, NekDouble> tmp;
1579  tmp = input[nlevels-1];
1580 
1581  for(int n = nlevels-1; n > 0; --n)
1582  {
1583  input[n] = input[n-1];
1584  }
1585 
1586  input[0] = tmp;
1587 }
void Nektar::ForcingMovingBody::SetDynEqCoeffMatrix ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields)
private

Definition at line 1467 of file ForcingMovingBody.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, m_CoeffMat_A, m_CoeffMat_B, m_lhom, Nektar::SolverUtils::Forcing::m_session, m_structdamp, m_structrho, and m_timestep.

Referenced by InitialiseCableModel().

1469 {
1470  int nplanes;
1471 
1472  bool homostrip;
1473  m_session->MatchSolverInfo("HomoStrip","True",homostrip,false);
1474 
1475  if(!homostrip)
1476  {
1477  nplanes = m_session->GetParameter("HomModesZ");
1478  }
1479  else
1480  {
1481  m_session->LoadParameter("Strip_Z", nplanes);
1482  }
1483 
1484  m_CoeffMat_A = Array<OneD, DNekMatSharedPtr>(nplanes);
1485  m_CoeffMat_B = Array<OneD, DNekMatSharedPtr>(nplanes);
1486 
1487  NekDouble tmp1, tmp2, tmp3;
1488  NekDouble tmp4, tmp5, tmp6, tmp7;
1489 
1490  // load the structural dynamic parameters from xml file
1491  NekDouble cabletension;
1492  NekDouble bendingstiff;
1493  NekDouble structstiff;
1494  m_session->LoadParameter("StructStiff", structstiff, 0.0);
1495  m_session->LoadParameter("CableTension", cabletension, 0.0);
1496  m_session->LoadParameter("BendingStiff", bendingstiff, 0.0);
1497 
1498  tmp1 = m_timestep * m_timestep;
1499  tmp2 = structstiff / m_structrho;
1500  tmp3 = m_structdamp / m_structrho;
1501  tmp4 = cabletension / m_structrho;
1502  tmp5 = bendingstiff / m_structrho;
1503 
1504  // solve the ODE in the wave space for cable motion to obtain disp, vel and
1505  // accel
1506 
1507  std::string supptype = m_session->GetSolverInfo("SupportType");
1508 
1509  for(int plane = 0; plane < nplanes; plane++)
1510  {
1511  int nel = 3;
1512  m_CoeffMat_A[plane]
1513  = MemoryManager<DNekMat>::AllocateSharedPtr(nel,nel);
1514  m_CoeffMat_B[plane]
1515  = MemoryManager<DNekMat>::AllocateSharedPtr(nel,nel);
1516 
1517  // Initialised to avoid compiler warnings.
1518  unsigned int K = 0;
1519  NekDouble beta = 0.0;
1520 
1521  if (boost::iequals(supptype, "Free-Free"))
1522  {
1523  K = plane/2;
1524  beta = 2.0 * M_PI/m_lhom;
1525  }
1526  else if(boost::iequals(supptype, "Pinned-Pinned"))
1527  {
1528  K = plane+1;
1529  beta = M_PI/m_lhom;
1530  }
1531  else
1532  {
1533  ASSERTL0(false,
1534  "Unrecognized support type for cable's motion");
1535  }
1536 
1537  tmp6 = beta * K;
1538  tmp6 = tmp6 * tmp6;
1539  tmp7 = tmp6 * tmp6;
1540  tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1541 
1542  (*m_CoeffMat_A[plane])(0,0) = tmp7;
1543  (*m_CoeffMat_A[plane])(0,1) = tmp3;
1544  (*m_CoeffMat_A[plane])(0,2) = 1.0;
1545  (*m_CoeffMat_A[plane])(1,0) = 0.0;
1546  (*m_CoeffMat_A[plane])(1,1) = 1.0;
1547  (*m_CoeffMat_A[plane])(1,2) =-m_timestep/2.0;
1548  (*m_CoeffMat_A[plane])(2,0) = 1.0;
1549  (*m_CoeffMat_A[plane])(2,1) = 0.0;
1550  (*m_CoeffMat_A[plane])(2,2) =-tmp1/4.0;
1551  (*m_CoeffMat_B[plane])(0,0) = 0.0;
1552  (*m_CoeffMat_B[plane])(0,1) = 0.0;
1553  (*m_CoeffMat_B[plane])(0,2) = 0.0;
1554  (*m_CoeffMat_B[plane])(1,0) = 0.0;
1555  (*m_CoeffMat_B[plane])(1,1) = 1.0;
1556  (*m_CoeffMat_B[plane])(1,2) = m_timestep/2.0;
1557  (*m_CoeffMat_B[plane])(2,0) = 1.0;
1558  (*m_CoeffMat_B[plane])(2,1) = m_timestep;
1559  (*m_CoeffMat_B[plane])(2,2) = tmp1/4.0;
1560 
1561  m_CoeffMat_A[plane]->Invert();
1562  (*m_CoeffMat_B[plane]) =
1563  (*m_CoeffMat_A[plane]) * (*m_CoeffMat_B[plane]);
1564  }
1565 }
Nektar::ForcingMovingBody::m_CoeffMat_B
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
Definition: ForcingMovingBody.h:157
Nektar::ForcingMovingBody::m_structrho
NekDouble m_structrho
mass of the cable per unit length
Definition: ForcingMovingBody.h:135
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ForcingMovingBody::m_structdamp
NekDouble m_structdamp
damping ratio of the cable
Definition: ForcingMovingBody.h:136
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::ForcingMovingBody::m_CoeffMat_A
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
Definition: ForcingMovingBody.h:155
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::ForcingMovingBody::m_timestep
NekDouble m_timestep
time step
Definition: ForcingMovingBody.h:139
Nektar::ForcingMovingBody::m_lhom
NekDouble m_lhom
length ratio of the cable
Definition: ForcingMovingBody.h:137
void Nektar::ForcingMovingBody::TensionedCableModel ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
Array< OneD, NekDouble > &  FcePhysinArray,
Array< OneD, NekDouble > &  MotPhysinArray 
)
private

Definition at line 580 of file ForcingMovingBody.cpp.

References ASSERTL0, m_CoeffMat_A, m_CoeffMat_B, m_FFT, m_np, Nektar::SolverUtils::Forcing::m_session, m_structrho, npts, and Vmath::Vcopy().

Referenced by EvaluateStructDynModel().

584 {
585  std::string supptype = m_session->GetSolverInfo("SupportType");
586 
587  bool homostrip;
588  m_session->MatchSolverInfo("HomoStrip","True",homostrip,false);
589  // homostrip == false indicates that a single fourier transformation is
590  // implemented for the motion of cable, so it is matched with that carried
591  // out in fluid fields; on the other hand, if homostrip == true, there is
592  // a mismatch between the structure and fluid fields in terms of fourier
593  // transformation, then the routines such as FFTFwdTrans/
594  // FFTBwdTrans can not be used dimectly for cable's solution.
595 
596  int npts;
597 
598  if(!homostrip) //full resolutions
599  {
600  npts = m_session->GetParameter("HomModesZ");
601  }
602  else
603  {
604  m_session->LoadParameter("Strip_Z", npts);
605  }
606 
607  Array<OneD, Array<OneD, NekDouble> > fft_i(4);
608  Array<OneD, Array<OneD, NekDouble> > fft_o(4);
609 
610  for(int i = 0; i < 4; i++)
611  {
612  fft_i[i] = Array<OneD, NekDouble>(npts, 0.0);
613  fft_o[i] = Array<OneD, NekDouble>(npts, 0.0);
614  }
615 
616  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
617  int colrank = vcomm->GetColumnComm()->GetRank();
618  int nproc = vcomm->GetColumnComm()->GetSize();
619 
620  if(!homostrip) //full resolutions
621  {
622  Array<OneD, NekDouble> tmp(4*m_np);
623 
624  for(int n = 0; n < m_np; n++)
625  {
626  tmp[n] = AeroForces[n];
627  for(int j = 0; j < 3; ++j)
628  {
629  tmp[(j+1)*m_np + n] =
630  CableMotions[j*m_np + n];
631  }
632  }
633  // Send to root process.
634  if(colrank == 0)
635  {
636  Vmath::Vcopy(m_np, AeroForces, 1, fft_i[0], 1);
637 
638  for (int var = 0; var < 3; var++)
639  {
640  Vmath::Vcopy(m_np, CableMotions+var*m_np, 1, fft_i[var+1], 1);
641  }
642 
643  for (int i = 1; i < nproc; ++i)
644  {
645  vcomm->GetColumnComm()->Recv(i, tmp);
646  for(int n = 0; n < m_np; n++)
647  {
648  for(int var = 0; var < 4; ++var)
649  {
650  fft_i[var][i*m_np + n] = tmp[var*m_np + n];
651  }
652  }
653  }
654  }
655  else
656  {
657  vcomm->GetColumnComm()->Send(0, tmp);
658  }
659  }
660  else //strip modeling
661  {
662  Array<OneD, NekDouble> tmp(4);
663 
664  tmp[0] = AeroForces[0];
665  for(int j = 0; j < 3; ++j)
666  {
667  tmp[j+1] = CableMotions[j*m_np];
668  }
669 
670  // Send to root process.
671  if(colrank == 0)
672  {
673  for(int j = 0 ; j < 4; ++j)
674  {
675  fft_i[j][0] = tmp[j];
676  }
677 
678  for(int i = 1; i < npts; ++i)
679  {
680  vcomm->GetColumnComm()->Recv(i, tmp);
681 
682  for(int j = 0 ; j < 4; ++j)
683  {
684  fft_i[j][i] = tmp[j];
685  }
686  }
687  }
688  else
689  {
690  for(int i = 1; i < npts; ++i)
691  {
692  if(colrank == i)
693  {
694  vcomm->GetColumnComm()->Send(0, tmp);
695  }
696  }
697  }
698  }
699 
700  if(colrank == 0)
701  {
702  // Implement Fourier transformation of the motion variables
703  if(boost::iequals(supptype, "Free-Free"))
704  {
705  for(int j = 0 ; j < 4; ++j)
706  {
707  m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
708  }
709  }
710  else if(boost::iequals(supptype, "Pinned-Pinned"))
711  {
712  //TODO:
713  int N = fft_i[0].num_elements();
714 
715  for(int var = 0; var < 4; var++)
716  {
717  for(int i = 0; i < N; i++)
718  {
719  fft_o[var][i] = 0;
720 
721  for(int k = 0; k < N; k++)
722  {
723  fft_o[var][i] +=
724  fft_i[var][k]*
725  sin(M_PI/(N)*(k+1/2)*(i+1));
726  }
727  }
728  }
729  }
730  else
731  {
732  ASSERTL0(false,
733  "Unrecognized support type for cable's motion");
734  }
735 
736  // solve the ODE in the wave space
737  for(int i = 0; i < npts; ++i)
738  {
739  int nrows = 3;
740 
741  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
742  tmp0 = Array<OneD, NekDouble> (3,0.0);
743  tmp1 = Array<OneD, NekDouble> (3,0.0);
744  tmp2 = Array<OneD, NekDouble> (3,0.0);
745 
746  for(int var = 0; var < 3; var++)
747  {
748  tmp0[var] = fft_o[var+1][i];
749  }
750 
751  tmp2[0] = fft_o[0][i];
752 
753  Blas::Dgemv('N', nrows, nrows, 1.0,
754  &(m_CoeffMat_B[i]->GetPtr())[0],
755  nrows, &tmp0[0], 1,
756  0.0, &tmp1[0], 1);
757  Blas::Dgemv('N', nrows, nrows, 1.0/m_structrho,
758  &(m_CoeffMat_A[i]->GetPtr())[0],
759  nrows, &tmp2[0], 1,
760  1.0, &tmp1[0], 1);
761 
762  for(int var = 0; var < 3; var++)
763  {
764  fft_i[var][i] = tmp1[var];
765  }
766  }
767 
768  // get physical coeffients via Backward fourier transformation of wave
769  // coefficients
770  if(boost::iequals(supptype, "Free-Free"))
771  {
772  for(int var = 0; var < 3; var++)
773  {
774  m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
775  }
776  }
777  else if(boost::iequals(supptype, "Pinned-Pinned"))
778  {
779  //TODO:
780  int N = fft_i[0].num_elements();
781 
782  for(int var = 0; var < 3; var++)
783  {
784  for(int i = 0; i < N; i++)
785  {
786  fft_o[var][i] = 0;
787 
788  for(int k = 0; k < N; k++)
789  {
790  fft_o[var][i] +=
791  fft_i[var][k]*
792  sin(M_PI/(N)*(k+1)*(i+1/2))*2/N;
793  }
794  }
795  }
796  }
797  else
798  {
799  ASSERTL0(false,
800  "Unrecognized support type for cable's motion");
801  }
802  }
803 
804  // send physical coeffients to all planes of each processor
805  if(!homostrip)//full resolutions
806  {
807  Array<OneD, NekDouble> tmp(3*m_np);
808 
809  if(colrank != 0)
810  {
811  vcomm->GetColumnComm()->Recv(0, tmp);
812  Vmath::Vcopy(3*m_np, tmp, 1,
813  CableMotions, 1);
814  }
815  else
816  {
817  for (int i = 1; i < nproc; ++i)
818  {
819  for(int j = 0; j < 3; j++)
820  {
821  for(int n = 0; n < m_np; n++)
822  {
823  tmp[j*m_np+n] = fft_o[j][i*m_np+n];
824  }
825  }
826  vcomm->GetColumnComm()->Send(i, tmp);
827  }
828 
829  for(int j = 0; j < 3; j++)
830  {
831  for(int n = 0; n < m_np; n++)
832  {
833  tmp[j*m_np+n] = fft_o[j][n];
834  }
835  Vmath::Vcopy(3*m_np, tmp, 1,
836  CableMotions, 1);
837  }
838  }
839  }
840  else //strip modeling
841  {
842  Array<OneD, NekDouble> tmp(4);
843 
844  if(colrank != 0)
845  {
846  for (int j = 1; j < nproc/npts; j++)
847  {
848  if(colrank == j*npts)
849  {
850  vcomm->GetColumnComm()->Recv(0, tmp);
851 
852  for(int plane = 0; plane < m_np; plane++)
853  {
854  for(int var = 0; var < 3; var++)
855  {
856  CableMotions[var*m_np+plane]= tmp[var];
857  }
858  }
859  }
860  }
861 
862  for(int i = 1; i < npts; i++)
863  {
864  for (int j = 0; j < nproc/npts; j++)
865  {
866  if(colrank == i+j*npts)
867  {
868  vcomm->GetColumnComm()->Recv(0, tmp);
869 
870  for(int plane = 0; plane < m_np; plane++)
871  {
872  for(int var = 0; var < 3; var++)
873  {
874  CableMotions[var*m_np+plane] = tmp[var];
875  }
876  }
877  }
878  }
879  }
880  }
881  else
882  {
883  for(int j = 0; j < 3; ++j)
884  {
885  tmp[j] = fft_o[j][0];
886  }
887 
888  for (int j = 1; j < nproc/npts; j++)
889  {
890  vcomm->GetColumnComm()->Send(j*npts, tmp);
891  }
892 
893  for(int plane = 0; plane < m_np; plane++)
894  {
895  for(int var = 0; var < 3; var++)
896  {
897  CableMotions[var*m_np+plane] = tmp[var];
898  }
899  }
900 
901  for(int i = 1; i < npts; ++i)
902  {
903  for(int j = 0; j < 3; ++j)
904  {
905  tmp[j] = fft_o[j][i];
906  }
907 
908  for (int j = 0; j < nproc/npts; j++)
909  {
910  vcomm->GetColumnComm()->Send(i+j*npts, tmp);
911  }
912  }
913  }
914  }
915 }
Nektar::ForcingMovingBody::m_CoeffMat_B
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
Definition: ForcingMovingBody.h:157
Nektar::ForcingMovingBody::m_structrho
NekDouble m_structrho
mass of the cable per unit length
Definition: ForcingMovingBody.h:135
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ForcingMovingBody::m_CoeffMat_A
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
Definition: ForcingMovingBody.h:155
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::ForcingMovingBody::m_FFT
LibUtilities::NektarFFTSharedPtr m_FFT
Definition: ForcingMovingBody.h:141
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::ForcingMovingBody::m_np
int m_np
number of planes per processors
Definition: ForcingMovingBody.h:132
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void Nektar::ForcingMovingBody::UpdateMotion ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const Array< OneD, Array< OneD, NekDouble > > &  fields,
NekDouble  time 
)
private

Definition at line 134 of file ForcingMovingBody.cpp.

References ASSERTL0, Nektar::MultiRegions::DirCartesianMap, Nektar::SolverUtils::Forcing::EvaluateFunction(), EvaluateStructDynModel(), m_Aeroforces, m_eta, m_funcName, m_IsFromFile, m_motion, m_MotionVars, m_MovBodyfilter, m_np, Nektar::SolverUtils::Forcing::m_session, m_zta, and Vmath::Vmul().

Referenced by v_Apply().

138 {
139  // Update the forces from the calculation of fluid field, which is
140  // implemented in the movingbody filter
141  m_MovBodyfilter->UpdateForce(m_session, pFields, m_Aeroforces, time);
142 
143  // for "free" type, the cable vibrates both in streamwise and crossflow
144  // dimections, for "constrained" type, the cable only vibrates in crossflow
145  // dimection, and for "forced" type, the calbe vibrates specifically along
146  // a given function or file
147  std::string vibtype = m_session->GetSolverInfo("VibrationType");
148  if(boost::iequals(vibtype, "Free") || boost::iequals(vibtype,"Constrained"))
149  {
150  // For free vibration case, displacements, velocities and acceleartions
151  // are obtained through solving structure dynamic model
152  EvaluateStructDynModel(pFields, time);
153  }
154  else if(boost::iequals(vibtype, "Forced"))
155  {
156  // For forced vibration case, load from specified file or function
157  int cnt = 0;
158  for(int j = 0; j < m_funcName.num_elements(); j++)
159  {
160  if(m_IsFromFile[cnt] && m_IsFromFile[cnt+1])
161  {
162  ASSERTL0(false, "Motion loading from file needs specific "
163  "implementation: Work in Progress!");
164  }
165  else
166  {
167  EvaluateFunction(pFields, m_session, m_motion[0], m_zta[j],
168  m_funcName[j], time);
169  EvaluateFunction(pFields, m_session, m_motion[1], m_eta[j],
170  m_funcName[j], time);
171  cnt = cnt + 2;
172  }
173  }
174 
175  for(int var = 0; var < 3; var++)
176  {
177  for(int plane = 0; plane < m_np; plane++)
178  {
179  int n = pFields[0]->GetPlane(plane)->GetTotPoints();
180  int offset = plane * n;
181  int xoffset = var * m_np+plane;
182  int yoffset = 3*m_np+xoffset;
183 
184  m_MotionVars[xoffset] = m_zta[var][offset];
185  m_MotionVars[yoffset] = m_eta[var][offset];
186  }
187  }
188  }
189  else
190  {
191  ASSERTL0(false,
192  "Unrecogized vibration type for cable's dynamic model");
193  }
194 
195  // Pass the variables of the cable's motion to the movingbody filter
196  m_MovBodyfilter->UpdateMotion(m_session, pFields, m_MotionVars, time);
197 
198  // Now we need to calcualte all the required z-derivatives
199  bool OriginalWaveSpace = pFields[0]->GetWaveSpace();
200  pFields[0]->SetWaveSpace(false);
201 
202  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
203  m_zta[0], m_zta[3]); //d(m_zta)/dz
204  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
205  m_zta[3], m_zta[4]); //dd(m_zta)/ddzz
206  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
207  m_zta[4], m_zta[5]); //ddd(m_zta)/dddzzz
208 
209  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
210  m_eta[0], m_eta[3]); //d(m_eta)/dz
211  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
212  m_eta[3], m_eta[4]); //dd(m_eta)/ddzz
213  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
214  m_eta[4], m_eta[5]); //ddd(m_eta)/dddzzz
215 
216  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
217  m_zta[1], m_zta[6]); //dd(m_zta)/ddtz
218  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
219  m_zta[6], m_zta[7]); //ddd(m_zta)/ddtzz
220 
221  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
222  m_eta[1], m_eta[6]); //dd(m_eta)/ddtz
223  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
224  m_eta[6], m_eta[7]); //ddd(m_eta)/ddtzz
225 
226  int NumPoints = pFields[0]->GetTotPoints();
227 
228  // (d(m_zta)/dz)(d(m_eta)/dz)
229  Vmath::Vmul(NumPoints, m_zta[3], 1, m_eta[3], 1, m_zta[8], 1);
230  // (d(m_eta)/dz)(d(m_zta)/dz) // not really needed
231  Vmath::Vmul(NumPoints, m_eta[3], 1, m_zta[3], 1, m_eta[8], 1);
232 
233  // (d(m_zta)/dz)^2
234  Vmath::Vmul(NumPoints, m_zta[3], 1, m_zta[3], 1, m_zta[9], 1);
235  // (d(m_eta)/dz)^2
236  Vmath::Vmul(NumPoints, m_eta[3], 1, m_eta[3], 1, m_eta[9], 1);
237 
238  pFields[0]->SetWaveSpace(OriginalWaveSpace);
239 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ForcingMovingBody::EvaluateStructDynModel
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, NekDouble time)
Definition: ForcingMovingBody.cpp:921
Nektar::ForcingMovingBody::m_Aeroforces
Array< OneD, NekDouble > m_Aeroforces
storage for the cable's force(x,y) variables
Definition: ForcingMovingBody.h:145
Nektar::SolverUtils::Forcing::EvaluateFunction
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, MultiRegions::ExpListSharedPtr > pFields, LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
Definition: Forcing.cpp:140
Nektar::ForcingMovingBody::m_motion
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
Definition: ForcingMovingBody.h:161
Nektar::ForcingMovingBody::m_MotionVars
Array< OneD, NekDouble > m_MotionVars
storage for the cable's motion(x,y) variables
Definition: ForcingMovingBody.h:147
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Nektar::ForcingMovingBody::m_IsFromFile
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
Definition: ForcingMovingBody.h:163
Nektar::ForcingMovingBody::m_eta
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
Definition: ForcingMovingBody.h:167
Nektar::ForcingMovingBody::m_MovBodyfilter
FilterMovingBodySharedPtr m_MovBodyfilter
Definition: ForcingMovingBody.h:143
Nektar::ForcingMovingBody::m_np
int m_np
number of planes per processors
Definition: ForcingMovingBody.h:132
Nektar::ForcingMovingBody::m_funcName
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
Definition: ForcingMovingBody.h:159
Nektar::ForcingMovingBody::m_zta
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
Definition: ForcingMovingBody.h:165
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::ForcingMovingBody::v_Apply ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const NekDouble &  time 
)
protectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 109 of file ForcingMovingBody.cpp.

References CalculateForcing(), m_forcing, MappingBndConditions(), UpdateMotion(), and Vmath::Vadd().

114 {
115  // Fill in m_zta and m_eta with all the required values
116  UpdateMotion(fields,inarray,time);
117 
118  //calcualte the m_forcing components Ax,Ay,Az and put them in m_forcing
119  CalculateForcing(fields);
120  // Apply m_forcing terms
121  for (int i = 0; i < 3; i++)
122  {
123  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1,
124  m_forcing[i], 1, outarray[i], 1);
125  }
126 
127  // Mapping boundary conditions
128  MappingBndConditions(fields, inarray, time);
129 }
Nektar::ForcingMovingBody::CalculateForcing
void CalculateForcing(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: ForcingMovingBody.cpp:245
Nektar::ForcingMovingBody::m_forcing
Array< OneD, Array< OneD, NekDouble > > m_forcing
Definition: ForcingMovingBody.h:169
Nektar::ForcingMovingBody::MappingBndConditions
void MappingBndConditions(const Array< OneD, MultiRegions::ExpListSharedPtr > &pfields, const Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
Definition: ForcingMovingBody.cpp:1028
Nektar::ForcingMovingBody::UpdateMotion
void UpdateMotion(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
Definition: ForcingMovingBody.cpp:134
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void Nektar::ForcingMovingBody::v_InitObject ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields,
const unsigned int &  pNumForcingFields,
const TiXmlElement *  pForce 
)
protectedvirtual

Implements Nektar::SolverUtils::Forcing.

Definition at line 53 of file ForcingMovingBody.cpp.

References ASSERTL0, CheckIsFromFile(), Nektar::MultiRegions::e3DH1D, InitialiseCableModel(), InitialiseFilter(), m_eta, m_forcing, Nektar::SolverUtils::Forcing::m_session, and m_zta.

57 {
58  // Just 3D homogenous 1D problems can use this techinque
59  ASSERTL0(pFields[0]->GetExpType()==MultiRegions::e3DH1D,
60  "Moving body implemented just for 3D Homogenous 1D expansions.");
61 
62  // At this point we know in the xml file where those quantities
63  // are declared (equation or file) - via a function name which is now
64  // stored in funcNameD etc. We need now to fill in with this info the
65  // m_zta and m_eta vectors (actuallythey are matrices) Array to control if
66  // the motion is determined by an equation or is from a file.(not Nektar++)
67  // check if we need to load a file or we have an equation
68  CheckIsFromFile(pForce);
69 
70  // Initialise movingbody filter
71  InitialiseFilter(m_session, pFields, pForce);
72 
73  // Initialise the cable model
74  InitialiseCableModel(m_session, pFields);
75 
76  m_zta = Array<OneD, Array< OneD, NekDouble> > (10);
77  m_eta = Array<OneD, Array< OneD, NekDouble> > (10);
78  // What are this bi-dimensional vectors ------------------------------------------
79  // m_zta[0] = m_zta | m_eta[0] = m_eta |
80  // m_zta[1] = d(m_zta)/dt | m_eta[1] = d(m_eta)/dt |
81  // m_zta[2] = dd(m_zta)/ddtt | m_eta[2] = dd(m_eta)/ddtt |
82  // m_zta[3] = d(m_zta)/dz | m_eta[3] = d(m_eta)/dz |
83  // m_zta[4] = dd(m_zta)/ddzz | m_eta[4] = dd(m_eta)/ddzz |
84  // m_zta[5] = ddd(m_zta)/dddzzz | m_eta[5] = ddd(m_eta)/dddzzz |
85  // m_zta[6] = dd(m_zta)/ddtz | m_eta[6] = dd(m_eta)/ddtz |
86  // m_zta[7] = ddd(m_zta)/dddtzz | m_eta[7] = ddd(m_eta)/dddtzz |
87  // m_zta[8] = (d(m_zta)/dz)(d(m_eta)/dz)| m_eta[8] = (d(m_eta)/dz)(d(m_zta)/dz) |
88  // m_zta[9] = (d(zm_zta)/dz)^2 | m_eta[9] = (d(m_eta)/dz)^2 |
89  //--------------------------------------------------------------------------------
90  int phystot = pFields[0]->GetTotPoints();
91  for(int i = 0; i < m_zta.num_elements(); i++)
92  {
93  m_zta[i] = Array<OneD, NekDouble>(phystot,0.0);
94  m_eta[i] = Array<OneD, NekDouble>(phystot,0.0);
95  }
96 
97  // m_forcing size (it must be 3: Ax, Ay and Az)
98  // total number of structural variables(Disp, Vel and Accel) must be 3
99  m_forcing = Array<OneD, Array<OneD, NekDouble> >(3);
100  // create the storage space for the body motion description
101 
102  for(int i = 0; i < 3; i++)
103  {
104  m_forcing[i] = Array<OneD, NekDouble> (phystot, 0.0);
105  }
106 
107 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ForcingMovingBody::CheckIsFromFile
void CheckIsFromFile(const TiXmlElement *pForce)
Definition: ForcingMovingBody.cpp:1628
Nektar::ForcingMovingBody::InitialiseFilter
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
Definition: ForcingMovingBody.cpp:1779
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::ForcingMovingBody::InitialiseCableModel
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Definition: ForcingMovingBody.cpp:1106
Nektar::ForcingMovingBody::m_eta
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
Definition: ForcingMovingBody.h:167
Nektar::ForcingMovingBody::m_forcing
Array< OneD, Array< OneD, NekDouble > > m_forcing
Definition: ForcingMovingBody.h:169
Nektar::ForcingMovingBody::m_zta
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
Definition: ForcingMovingBody.h:165

Friends And Related Function Documentation

friend class MemoryManager< ForcingMovingBody >
friend

Definition at line 53 of file ForcingMovingBody.h.

Member Data Documentation

std::string Nektar::ForcingMovingBody::className
static
Initial value:
= SolverUtils::GetForcingFactory().
RegisterCreatorFunction("MovingBody",
ForcingMovingBody::create,
"Moving Body Forcing")

Name of the class.

Definition at line 70 of file ForcingMovingBody.h.

Array<OneD, NekDouble> Nektar::ForcingMovingBody::m_Aeroforces
private

storage for the cable's force(x,y) variables

Definition at line 145 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), and UpdateMotion().

Array<OneD, DNekMatSharedPtr> Nektar::ForcingMovingBody::m_CoeffMat_A
private

matrices in Newmart-beta method

Definition at line 155 of file ForcingMovingBody.h.

Referenced by SetDynEqCoeffMatrix(), and TensionedCableModel().

Array<OneD, DNekMatSharedPtr> Nektar::ForcingMovingBody::m_CoeffMat_B
private

matrices in Newmart-beta method

Definition at line 157 of file ForcingMovingBody.h.

Referenced by SetDynEqCoeffMatrix(), and TensionedCableModel().

Array<OneD, Array< OneD, NekDouble> > Nektar::ForcingMovingBody::m_eta
private

Store the derivatives of motion variables in y-direction.

Definition at line 167 of file ForcingMovingBody.h.

Referenced by CalculateForcing(), EvaluateStructDynModel(), UpdateMotion(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::ForcingMovingBody::m_fA
private

fictitious acceleration storage

Definition at line 153 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

LibUtilities::NektarFFTSharedPtr Nektar::ForcingMovingBody::m_FFT
private

Definition at line 141 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and TensionedCableModel().

Array<OneD, Array< OneD, NekDouble> > Nektar::ForcingMovingBody::m_forcing
private

Definition at line 169 of file ForcingMovingBody.h.

Referenced by CalculateForcing(), v_Apply(), and v_InitObject().

Array<OneD, std::string> Nektar::ForcingMovingBody::m_funcName
private

[0] is displacements, [1] is velocities, [2] is accelerations

Definition at line 159 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), InitialiseCableModel(), and UpdateMotion().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::ForcingMovingBody::m_fV
private

fictitious velocity storage

Definition at line 151 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

Array<OneD, bool> Nektar::ForcingMovingBody::m_IsFromFile
private

do determine if the the body motion come from an extern file

Definition at line 163 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), InitialiseCableModel(), and UpdateMotion().

NekDouble Nektar::ForcingMovingBody::m_kinvis
private

fluid viscous

Definition at line 138 of file ForcingMovingBody.h.

Referenced by CalculateForcing(), and InitialiseCableModel().

NekDouble Nektar::ForcingMovingBody::m_lhom
private

length ratio of the cable

Definition at line 137 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), MappingBndConditions(), and SetDynEqCoeffMatrix().

Array<OneD, std::string> Nektar::ForcingMovingBody::m_motion
private

motion direction: [0] is 'x' and [1] is 'y'

Definition at line 161 of file ForcingMovingBody.h.

Referenced by CheckIsFromFile(), EvaluateStructDynModel(), InitialiseCableModel(), MappingBndConditions(), and UpdateMotion().

Array<OneD, NekDouble> Nektar::ForcingMovingBody::m_MotionVars
private

storage for the cable's motion(x,y) variables

Definition at line 147 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), MappingBndConditions(), and UpdateMotion().

FilterMovingBodySharedPtr Nektar::ForcingMovingBody::m_MovBodyfilter
private

Definition at line 143 of file ForcingMovingBody.h.

Referenced by InitialiseFilter(), and UpdateMotion().

int Nektar::ForcingMovingBody::m_movingBodyCalls
private

number of times the movbody have been called

Definition at line 131 of file ForcingMovingBody.h.

Referenced by EvaluateAccelaration(), EvaluateStructDynModel(), and InitialiseCableModel().

int Nektar::ForcingMovingBody::m_np
private

number of planes per processors

Definition at line 132 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), InitialiseCableModel(), MappingBndConditions(), TensionedCableModel(), and UpdateMotion().

NekDouble Nektar::ForcingMovingBody::m_structdamp
private

damping ratio of the cable

Definition at line 136 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), and SetDynEqCoeffMatrix().

NekDouble Nektar::ForcingMovingBody::m_structrho
private

mass of the cable per unit length

Definition at line 135 of file ForcingMovingBody.h.

Referenced by InitialiseCableModel(), SetDynEqCoeffMatrix(), and TensionedCableModel().

NekDouble Nektar::ForcingMovingBody::m_timestep
private

time step

Definition at line 139 of file ForcingMovingBody.h.

Referenced by EvaluateAccelaration(), InitialiseCableModel(), and SetDynEqCoeffMatrix().

int Nektar::ForcingMovingBody::m_vdim
private

vibration dimension

Definition at line 133 of file ForcingMovingBody.h.

Referenced by EvaluateStructDynModel(), and InitialiseCableModel().

Array<OneD, Array<OneD, NekDouble> > Nektar::ForcingMovingBody::m_W
private

srorage for the velocity in z-direction

Definition at line 149 of file ForcingMovingBody.h.

Referenced by EvaluateAccelaration(), and InitialiseCableModel().

Array<OneD, Array< OneD, NekDouble> > Nektar::ForcingMovingBody::m_zta
private

Store the derivatives of motion variables in x-direction.

Definition at line 165 of file ForcingMovingBody.h.

Referenced by CalculateForcing(), EvaluateStructDynModel(), UpdateMotion(), and v_InitObject().

Generated on Tue Dec 8 2015 23:45:34 for Nektar++ by   doxygen 1.8.9.1