Nektar++
ForcingMovingBody.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ForcingMovingBody.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Moving Body m_forcing (movement of a body in a domain is achieved
33 // via a m_forcing term which is the results of a coordinate system motion)
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <IncNavierStokesSolver/Forcing/ForcingMovingBody.h>
38 #include <MultiRegions/ExpList.h>
39 
40 namespace Nektar
41 {
42 std::string ForcingMovingBody::className = SolverUtils::GetForcingFactory().
43  RegisterCreatorFunction("MovingBody",
44  ForcingMovingBody::create,
45  "Moving Body Forcing");
46 
47 ForcingMovingBody::ForcingMovingBody(
48  const LibUtilities::SessionReaderSharedPtr& pSession)
49  : Forcing(pSession)
50 {
51 }
52 
53 void ForcingMovingBody::v_InitObject(
54  const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
55  const unsigned int& pNumForcingFields,
56  const TiXmlElement* pForce)
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 }
108 
109 void ForcingMovingBody::v_Apply(
110  const Array<OneD, MultiRegions::ExpListSharedPtr>& fields,
111  const Array<OneD, Array<OneD, NekDouble> >& inarray,
112  Array<OneD, Array<OneD, NekDouble> >& outarray,
113  const NekDouble& time)
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 }
130 
131 /**
132  *
133  */
134 void ForcingMovingBody::UpdateMotion(
135  const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
136  const Array<OneD, Array<OneD, NekDouble> > & fields,
137  NekDouble time)
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 }
240 
241 
242 /**
243  *
244  */
245 void ForcingMovingBody::CalculateForcing(
246  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
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 }
575 
576 
577 /**
578  *
579  */
580 void ForcingMovingBody::TensionedCableModel(
581  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
582  Array<OneD, NekDouble> &AeroForces,
583  Array<OneD, NekDouble> &CableMotions)
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 }
916 
917 
918 /**
919  *
920  */
921 void ForcingMovingBody::EvaluateStructDynModel(
922  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
923  NekDouble time)
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 }
1024 
1025 /**
1026  *
1027  */
1028 void ForcingMovingBody::MappingBndConditions(
1029  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
1030  const Array<OneD, Array<OneD, NekDouble> > &fields,
1031  NekDouble time)
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 }
1101 
1102 
1103 /**
1104  *
1105  */
1106 void ForcingMovingBody::InitialiseCableModel(
1107  const LibUtilities::SessionReaderSharedPtr& pSession,
1108  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
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 }
1462 
1463 
1464 /**
1465  *
1466  */
1467 void ForcingMovingBody::SetDynEqCoeffMatrix(
1468  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
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 }
1566 
1567 
1568 /**
1569  * Function to roll time-level storages to the next step layout.
1570  * The stored data associated with the oldest time-level
1571  * (not required anymore) are moved to the top, where they will
1572  * be overwritten as the solution process progresses.
1573  */
1574 void ForcingMovingBody::RollOver(Array<OneD, Array<OneD, NekDouble> > &input)
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 }
1588 
1589 
1590 /**
1591  *
1592  */
1593 void ForcingMovingBody::EvaluateAccelaration(
1594  const Array<OneD, NekDouble> &input,
1595  Array<OneD, NekDouble> &output,
1596  int npoints)
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 }
1623 
1624 
1625 /**
1626  *
1627  */
1628 void ForcingMovingBody::CheckIsFromFile(const TiXmlElement* pForce)
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 }
1775 
1776 /**
1777  *
1778  */
1779 void ForcingMovingBody::InitialiseFilter(
1780  const LibUtilities::SessionReaderSharedPtr& pSession,
1781  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
1782  const TiXmlElement* pForce)
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 }
1815 
1816 }
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:135
Nektar::ForcingMovingBody::CheckIsFromFile
void CheckIsFromFile(const TiXmlElement *pForce)
Definition: ForcingMovingBody.cpp:1628
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::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::ForcingMovingBody::EvaluateStructDynModel
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, NekDouble time)
Definition: ForcingMovingBody.cpp:921
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
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_CoeffMat_A
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
Definition: ForcingMovingBody.h:155
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::SolverUtils::GetForcingFactory
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:42
Nektar::ForcingMovingBody::m_Aeroforces
Array< OneD, NekDouble > m_Aeroforces
storage for the cable's force(x,y) variables
Definition: ForcingMovingBody.h:145
Nektar::ForcingMovingBody::v_InitObject
virtual void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Definition: ForcingMovingBody.cpp:53
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::CalculateForcing
void CalculateForcing(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: ForcingMovingBody.cpp:245
Nektar::LibUtilities::GetNektarFFTFactory
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
Nektar::ForcingMovingBody::InitialiseFilter
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
Definition: ForcingMovingBody.cpp:1779
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
Nektar::ForcingMovingBody::ForcingMovingBody
ForcingMovingBody(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: ForcingMovingBody.cpp:47
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::ForcingMovingBody::EvaluateAccelaration
void EvaluateAccelaration(const Array< OneD, NekDouble > &input, Array< OneD, NekDouble > &output, int npoints)
Definition: ForcingMovingBody.cpp:1593
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
ForcingMovingBody.h
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
ExpList.h
Nektar::ForcingMovingBody::v_Apply
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)
Definition: ForcingMovingBody.cpp:109
Nektar::LibUtilities::Equation::Evaluate
NekDouble Evaluate() const
Definition: Equation.h:102
Nektar::ForcingMovingBody::TensionedCableModel
void TensionedCableModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
Definition: ForcingMovingBody.cpp:580
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::Forcing::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition: NektarUnivConsts.hpp:47
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
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::InitialiseCableModel
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Definition: ForcingMovingBody.cpp:1106
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_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::className
static std::string className
Name of the class.
Definition: ForcingMovingBody.h:70
Nektar::ForcingMovingBody::m_np
int m_np
number of planes per processors
Definition: ForcingMovingBody.h:132
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
Nektar::ForcingMovingBody::m_funcName
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
Definition: ForcingMovingBody.h:159
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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
Nektar::ForcingMovingBody::create
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.
Definition: ForcingMovingBody.h:56
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1016
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
Generated on Tue Dec 8 2015 23:41:28 for Nektar++ by   doxygen 1.8.9.1