Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
38 #include <MultiRegions/ExpList.h>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 std::string ForcingMovingBody::className = SolverUtils::GetForcingFactory().
45  RegisterCreatorFunction("MovingBody",
46  ForcingMovingBody::create,
47  "Moving Body Forcing");
48 
49 ForcingMovingBody::ForcingMovingBody(
51  : Forcing(pSession)
52 {
53 }
54 
57  const unsigned int& pNumForcingFields,
58  const TiXmlElement* pForce)
59 {
60  // Just 3D homogenous 1D problems can use this techinque
61  ASSERTL0(pFields[0]->GetExpType()==MultiRegions::e3DH1D,
62  "Moving body implemented just for 3D Homogenous 1D expansions.");
63 
64  // At this point we know in the xml file where those quantities
65  // are declared (equation or file) - via a function name which is now
66  // stored in funcNameD etc. We need now to fill in with this info the
67  // m_zta and m_eta vectors (actuallythey are matrices) Array to control if
68  // the motion is determined by an equation or is from a file.(not Nektar++)
69  // check if we need to load a file or we have an equation
70  CheckIsFromFile(pForce);
71 
72  // Initialise movingbody filter
73  InitialiseFilter(m_session, pFields, pForce);
74 
75  // Initialise the cable model
77 
78  // Load mapping
80  m_mapping->SetTimeDependent( true );
81 
82  if(m_vdim > 0)
83  {
84  m_mapping->SetFromFunction( false );
85  }
86 
89  // What are this bi-dimensional vectors ------------------------------------------
90  // m_zta[0] = m_zta | m_eta[0] = m_eta |
91  // m_zta[1] = d(m_zta)/dt | m_eta[1] = d(m_eta)/dt |
92  // m_zta[2] = dd(m_zta)/ddtt | m_eta[2] = dd(m_eta)/ddtt |
93  //--------------------------------------------------------------------------------
94  int phystot = pFields[0]->GetTotPoints();
95  for(int i = 0; i < m_zta.num_elements(); i++)
96  {
97  m_zta[i] = Array<OneD, NekDouble>(phystot,0.0);
98  m_eta[i] = Array<OneD, NekDouble>(phystot,0.0);
99  }
100 }
101 
104  const Array<OneD, Array<OneD, NekDouble> >& inarray,
105  Array<OneD, Array<OneD, NekDouble> >& outarray,
106  const NekDouble& time)
107 {
108  // Update the forces from the calculation of fluid field, which is
109  // implemented in the movingbody filter
110  Array<OneD, NekDouble> Hydroforces (2*m_np,0.0);
111  m_MovBodyfilter->UpdateForce(m_session, pFields, Hydroforces, time);
112 
113  // for "free" type (m_vdim = 2), the cable vibrates both in streamwise and crossflow
114  // dimections, for "constrained" type (m_vdim = 1), the cable only vibrates in crossflow
115  // dimection, and for "forced" type (m_vdim = 0), the calbe vibrates specifically along
116  // a given function or file
117  if(m_vdim == 1 || m_vdim == 2)
118  {
119  // For free vibration case, displacements, velocities and acceleartions
120  // are obtained through solving structure dynamic model
121  EvaluateStructDynModel(pFields, Hydroforces, time);
122 
123  // Convert result to format required by mapping
124  int physTot = pFields[0]->GetTotPoints();
127  for(int i =0; i<3; i++)
128  {
129  coords[i] = Array< OneD, NekDouble> (physTot, 0.0);
130  coordsVel[i] = Array< OneD, NekDouble> (physTot, 0.0);
131  }
132  // Get original coordinates
133  pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
134 
135  // Add displacement to coordinates
136  Vmath::Vadd(physTot, coords[0], 1, m_zta[0], 1, coords[0], 1);
137  Vmath::Vadd(physTot, coords[1], 1, m_eta[0], 1, coords[1], 1);
138  // Copy velocities
139  Vmath::Vcopy(physTot, m_zta[1], 1, coordsVel[0], 1);
140  Vmath::Vcopy(physTot, m_eta[1], 1, coordsVel[1], 1);
141 
142  // Update mapping
143  m_mapping->UpdateMapping(time, coords, coordsVel);
144  }
145  else if(m_vdim == 0)
146  {
147  // For forced vibration case, load from specified file or function
148  int cnt = 0;
149  for(int j = 0; j < m_funcName.num_elements(); j++)
150  {
151  if(m_IsFromFile[cnt] && m_IsFromFile[cnt+1])
152  {
153  ASSERTL0(false, "Motion loading from file needs specific "
154  "implementation: Work in Progress!");
155  }
156  else
157  {
158  EvaluateFunction(pFields, m_session, m_motion[0], m_zta[j],
159  m_funcName[j], time);
160  EvaluateFunction(pFields, m_session, m_motion[1], m_eta[j],
161  m_funcName[j], time);
162  cnt = cnt + 2;
163  }
164  }
165 
166  // Update mapping
167  m_mapping->UpdateMapping(time);
168 
169  // Convert result from mapping
170  int physTot = pFields[0]->GetTotPoints();
173  for(int i =0; i<3; i++)
174  {
175  coords[i] = Array< OneD, NekDouble> (physTot, 0.0);
176  coordsVel[i] = Array< OneD, NekDouble> (physTot, 0.0);
177  }
178  // Get original coordinates
179  pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
180 
181  // Get Coordinates and coord velocity from mapping
182  m_mapping->GetCartesianCoordinates(m_zta[0], m_eta[0], coords[2]);
183  m_mapping->GetCoordVelocity(coordsVel);
184 
185  // Calculate displacement
186  Vmath::Vsub(physTot, m_zta[0], 1, coords[0], 1, m_zta[0], 1);
187  Vmath::Vsub(physTot, m_eta[0], 1, coords[1], 1, m_eta[0], 1);
188 
189  Vmath::Vcopy(physTot, coordsVel[0], 1, m_zta[1], 1);
190  Vmath::Vcopy(physTot, coordsVel[1], 1, m_eta[1], 1);
191 
192  for(int var = 0; var < 3; var++)
193  {
194  for(int plane = 0; plane < m_np; plane++)
195  {
196  int n = pFields[0]->GetPlane(plane)->GetTotPoints();
197  int offset = plane * n;
198  int Offset = var * m_np+plane;
199 
200  m_MotionVars[0][Offset] = m_zta[var][offset];
201  m_MotionVars[1][Offset] = m_eta[var][offset];
202  }
203  }
204  }
205  else
206  {
207  ASSERTL0(false,
208  "Unrecogized vibration type for cable's dynamic model");
209  }
210 
211  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
212  int colrank = vcomm->GetColumnComm()->GetRank();
213  // Pass the variables of the cable's motion to the movingbody filter
214  if(colrank == 0)
215  {
216  int n = m_MotionVars[0].num_elements();
217  Array<OneD, NekDouble> tmpArray(2*n),tmp(n);
218  Vmath::Vcopy(n,m_MotionVars[0],1,tmpArray,1);
219  Vmath::Vcopy(n,m_MotionVars[1],1,tmp=tmpArray+n,1);
220  m_MovBodyfilter->UpdateMotion(m_session, pFields, tmpArray, time);
221  }
222 }
223 
224 
225 
226 /**
227  *
228  */
231  Array<OneD, NekDouble> &Hydroforces,
232  NekDouble time)
233 {
234  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
235  int colrank = vcomm->GetColumnComm()->GetRank();
236  int nproc = vcomm->GetColumnComm()->GetSize();
237 
238  bool homostrip;
239  m_session->MatchSolverInfo("HomoStrip","True",homostrip,false);
240 
241  //number of structral modes and number of strips
242  int npts, nstrips;
243  if(!homostrip) //full resolutions
244  {
245  npts = m_session->GetParameter("HomModesZ");
246  }
247  else
248  {
249  m_session->LoadParameter("HomStructModesZ", npts);
250  m_session->LoadParameter("Strip_Z", nstrips);
251 
252  //ASSERTL0(nstrips < = npts,
253  // "the number of struct. modes must be larger than that of the strips.");
254  }
255 
256  //the hydrodynamic forces
258  //forces in x-direction
259  fces[0] = Array <OneD, NekDouble> (npts,0.0);
260  //forces in y-direction
261  fces[1] = Array <OneD, NekDouble> (npts,0.0);
262 
263  //fill the force vectors
264  if(colrank == 0)
265  {
266  if(!homostrip) //full resolutions
267  {
268  Vmath::Vcopy(m_np, Hydroforces, 1, fces[0], 1);
269  Vmath::Vcopy(m_np, Hydroforces+m_np, 1, fces[1], 1);
270  }
271  else //strip modelling
272  {
273  fces[0][0] = Hydroforces[0];
274  fces[1][0] = Hydroforces[m_np];
275  }
276 
277  if(!homostrip) //full resolutions
278  {
280  for (int i = 1; i < nproc; ++i)
281  {
282  vcomm->GetColumnComm()->Recv(i, tmp);
283  for(int n = 0; n < m_np; n++)
284  {
285  for(int j = 0; j < 2; ++j)
286  {
287  fces[j][i*m_np + n] = tmp[j*m_np + n];
288  }
289  }
290  }
291  }
292  else //strip modelling
293  //if the body is submerged partly, the fces are filled partly
294  //by the flow induced forces
295  {
296  Array<OneD, NekDouble> tmp(2);
297  for(int i = 1; i < nstrips; ++i)
298  {
299  vcomm->GetColumnComm()->Recv(i, tmp);
300 
301  for(int j = 0 ; j < 2; ++j)
302  {
303  fces[j][i] = tmp[j];
304  }
305  }
306  }
307  }
308  else
309  {
310  if(!homostrip) //full resolutions
311  {
312  vcomm->GetColumnComm()->Send(0, Hydroforces);
313  }
314  else //strip modelling
315  {
316  for(int i = 1; i < nstrips; ++i)
317  {
318  if(colrank == i)
319  {
320  Array<OneD, NekDouble> tmp(2);
321  tmp[0] = Hydroforces[0];
322  tmp[1] = Hydroforces[m_np];
323  vcomm->GetColumnComm()->Send(0, tmp);
324  }
325  }
326  }
327  }
328 
329  if(colrank == 0)
330  {
331  // Fictitious mass method used to stablize the explicit coupling at
332  // relatively lower mass ratio
333  bool fictmass;
334  m_session->MatchSolverInfo("FictitiousMassMethod", "True",
335  fictmass, false);
336  if(fictmass)
337  {
338  NekDouble fictrho, fictdamp;
339  m_session->LoadParameter("FictMass", fictrho);
340  m_session->LoadParameter("FictDamp", fictdamp);
341 
342  static NekDouble Betaq_Coeffs[2][2] =
343  {{1.0, 0.0},{2.0, -1.0}};
344 
345  // only consider second order approximation for fictitious variables
346  int intOrder= 2;
347  int nint = min(m_movingBodyCalls+1,intOrder);
348  int nlevels = m_fV[0].num_elements();
349 
350  for(int i = 0; i < m_motion.num_elements(); ++i)
351  {
352  RollOver(m_fV[i]);
353  RollOver(m_fA[i]);
354 
355  int Voffset = npts;
356  int Aoffset = 2*npts;
357 
358  Vmath::Vcopy(npts, m_MotionVars[i]+Voffset, 1, m_fV[i][0], 1);
359  Vmath::Vcopy(npts, m_MotionVars[i]+Aoffset, 1, m_fA[i][0], 1);
360 
361  // Extrapolate to n+1
362  Vmath::Smul(npts,
363  Betaq_Coeffs[nint-1][nint-1],
364  m_fV[i][nint-1], 1,
365  m_fV[i][nlevels-1], 1);
366  Vmath::Smul(npts,
367  Betaq_Coeffs[nint-1][nint-1],
368  m_fA[i][nint-1], 1,
369  m_fA[i][nlevels-1], 1);
370 
371  for(int n = 0; n < nint-1; ++n)
372  {
373  Vmath::Svtvp(npts,
374  Betaq_Coeffs[nint-1][n],
375  m_fV[i][n],1,m_fV[i][nlevels-1],1,
376  m_fV[i][nlevels-1],1);
377  Vmath::Svtvp(npts,
378  Betaq_Coeffs[nint-1][n],
379  m_fA[i][n],1,m_fA[i][nlevels-1],1,
380  m_fA[i][nlevels-1],1);
381  }
382 
383  // Add the fictitious forces on the RHS of the equation
384  Vmath::Svtvp(npts, fictdamp,m_fV[i][nlevels-1],1,
385  fces[i],1,fces[i],1);
386  Vmath::Svtvp(npts, fictrho, m_fA[i][nlevels-1],1,
387  fces[i],1,fces[i],1);
388  }
389  }
390  }
391  //structural solver is implemented on the root process
392  if(colrank == 0)
393  {
394  // Tensioned cable model is evaluated in wave space
395  for(int n = 0, cn = 1; n < m_vdim; n++, cn--)
396  {
397  Newmark_betaSolver(pFields, fces[cn], m_MotionVars[cn]);
398  }
399  }
400 
401  Array<OneD, NekDouble> Motvars(2*2*m_np);
402  // send physical coeffients to all planes of each processor
403  if(!homostrip)//full resolutions
404  {
405  Array<OneD, NekDouble> tmp(2*2*m_np);
406 
407  if(colrank != 0)
408  {
409  vcomm->GetColumnComm()->Recv(0, tmp);
410  Vmath::Vcopy(2*2*m_np, tmp, 1, Motvars, 1);
411  }
412  else
413  {
414  for (int i = 1; i < nproc; ++i)
415  {
416  for(int j = 0; j < 2; j++) //moving dimensions
417  {
418  for(int k = 0; k < 2; k++) //disp. and vel.
419  {
420  for (int n = 0; n < m_np; n++)
421  {
422  tmp[j*2*m_np+k*m_np+n] = m_MotionVars[j][k*npts+i*m_np+n];
423  }
424  }
425  }
426  vcomm->GetColumnComm()->Send(i, tmp);
427  }
428 
429  for(int j = 0; j < 2; j++)
430  {
431  for(int k = 0; k < 2; k++)
432  {
433  for(int n = 0; n < m_np; n++)
434  {
435  tmp[j*2*m_np+k*m_np+n] = m_MotionVars[j][k*npts+n];
436  }
437  Vmath::Vcopy(2*2*m_np, tmp, 1, Motvars, 1);
438  }
439  }
440  }
441  }
442  else //strip modelling
443  {
444  Array<OneD, NekDouble> tmp(2*2);
445 
446  if(colrank != 0)
447  {
448  for (int j = 1; j < nproc/nstrips; j++)
449  {
450  if(colrank == j*nstrips)
451  {
452  vcomm->GetColumnComm()->Recv(0, tmp);
453 
454  for(int plane = 0; plane < m_np; plane++)
455  {
456  for(int var = 0; var < 2; var++)
457  {
458  for(int k = 0; k < 2; k++)
459  {
460  Motvars[var*2*m_np+k*m_np+plane]= tmp[var*2+k];
461  }
462  }
463  }
464  }
465  }
466 
467  for(int i = 1; i < nstrips; i++)
468  {
469  for (int j = 0; j < nproc/nstrips; j++)
470  {
471  if(colrank == i+j*nstrips)
472  {
473  vcomm->GetColumnComm()->Recv(0, tmp);
474 
475  for(int plane = 0; plane < m_np; plane++)
476  {
477  for(int var = 0; var < 2; var++)
478  {
479  for(int k = 0; k < 2; k++)
480  {
481  Motvars[var*2*m_np+k*m_np+plane] = tmp[var*2+k];
482  }
483  }
484  }
485  }
486  }
487  }
488  }
489  else
490  {
491  for(int j = 0; j < 2; ++j)
492  {
493  for(int k = 0; k < 2; ++k)
494  {
495  tmp[j*2+k] = m_MotionVars[j][k*npts];
496  }
497  }
498 
499  for (int j = 1; j < nproc/nstrips; j++)
500  {
501  vcomm->GetColumnComm()->Send(j*nstrips, tmp);
502  }
503 
504  for(int plane = 0; plane < m_np; plane++)
505  {
506  for(int j = 0; j < 2; j++)
507  {
508  for(int k = 0; k < 2; ++k)
509  {
510  Motvars[j*2*m_np+k*m_np+plane] = m_MotionVars[j][k*npts];
511  }
512  }
513  }
514 
515  for(int i = 1; i < nstrips; ++i)
516  {
517  for(int j = 0; j < 2; ++j)
518  {
519  for(int k = 0; k < 2; ++k)
520  {
521  tmp[j*2+k] = m_MotionVars[j][i+k*npts];
522  }
523  }
524 
525  for (int j = 0; j < nproc/nstrips; j++)
526  {
527  vcomm->GetColumnComm()->Send(i+j*nstrips, tmp);
528  }
529  }
530  }
531  }
532 
533  // Set the m_forcing term based on the motion of the cable
534  for(int var = 0; var < 2; var++)
535  {
536  for(int plane = 0; plane < m_np; plane++)
537  {
538  int n = pFields[0]->GetPlane(plane)->GetTotPoints();
539 
541 
542  int offset = plane * n;
543  int xoffset = var * m_np+plane;
544  int yoffset = 2*m_np + xoffset;
545 
546  Vmath::Fill(n, Motvars[xoffset], tmp = m_zta[var] + offset, 1);
547  Vmath::Fill(n, Motvars[yoffset], tmp = m_eta[var] + offset, 1);
548  }
549  }
550 }
551 
552 
553 /**
554  *
555  */
558  Array<OneD, NekDouble> &HydroForces,
559  Array<OneD, NekDouble> &BodyMotions)
560 {
561  std::string supptype = m_session->GetSolverInfo("SupportType");
562 
563  int npts = HydroForces.num_elements();
564 
567 
568  for(int i = 0; i < 4; i++)
569  {
570  fft_i[i] = Array<OneD, NekDouble>(npts, 0.0);
571  fft_o[i] = Array<OneD, NekDouble>(npts, 0.0);
572  }
573 
574  Vmath::Vcopy(npts, HydroForces, 1, fft_i[0], 1);
575  for(int i = 0; i < 3; i++)
576  {
577  Vmath::Vcopy(npts, BodyMotions+i*npts, 1, fft_i[i+1], 1);
578  }
579 
580  // Implement Fourier transformation of the motion variables
581  if(boost::iequals(supptype, "Free-Free"))
582  {
583  for(int j = 0 ; j < 4; ++j)
584  {
585  m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
586  }
587  }
588  else if(boost::iequals(supptype, "Pinned-Pinned"))
589  {
590  //TODO:
591  int N = fft_i[0].num_elements();
592 
593  for(int var = 0; var < 4; var++)
594  {
595  for(int i = 0; i < N; i++)
596  {
597  fft_o[var][i] = 0;
598 
599  for(int k = 0; k < N; k++)
600  {
601  fft_o[var][i] +=
602  fft_i[var][k]*
603  sin(M_PI/(N)*(k+1/2)*(i+1));
604  }
605  }
606  }
607  }
608  else
609  {
610  ASSERTL0(false,
611  "Unrecognized support type for cable's motion");
612  }
613 
614  // solve the ODE in the wave space
615  for(int i = 0; i < npts; ++i)
616  {
617  int nrows = 3;
618 
619  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
620  tmp0 = Array<OneD, NekDouble> (3,0.0);
621  tmp1 = Array<OneD, NekDouble> (3,0.0);
622  tmp2 = Array<OneD, NekDouble> (3,0.0);
623 
624  for(int var = 0; var < 3; var++)
625  {
626  tmp0[var] = fft_o[var+1][i];
627  }
628 
629  tmp2[0] = fft_o[0][i];
630 
631  Blas::Dgemv('N', nrows, nrows, 1.0,
632  &(m_CoeffMat_B[i]->GetPtr())[0],
633  nrows, &tmp0[0], 1,
634  0.0, &tmp1[0], 1);
635  Blas::Dgemv('N', nrows, nrows, 1.0/m_structrho,
636  &(m_CoeffMat_A[i]->GetPtr())[0],
637  nrows, &tmp2[0], 1,
638  1.0, &tmp1[0], 1);
639 
640  for(int var = 0; var < 3; var++)
641  {
642  fft_i[var][i] = tmp1[var];
643  }
644  }
645 
646  // get physical coeffients via Backward fourier transformation of wave
647  // coefficients
648  if(boost::iequals(supptype, "Free-Free"))
649  {
650  for(int var = 0; var < 3; var++)
651  {
652  m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
653  }
654  }
655  else if(boost::iequals(supptype, "Pinned-Pinned"))
656  {
657  //TODO:
658  int N = fft_i[0].num_elements();
659 
660  for(int var = 0; var < 3; var++)
661  {
662  for(int i = 0; i < N; i++)
663  {
664  fft_o[var][i] = 0;
665 
666  for(int k = 0; k < N; k++)
667  {
668  fft_o[var][i] +=
669  fft_i[var][k]*
670  sin(M_PI/(N)*(k+1)*(i+1/2))*2/N;
671  }
672  }
673  }
674  }
675  else
676  {
677  ASSERTL0(false,
678  "Unrecognized support type for cable's motion");
679  }
680 
681 
682  for(int i = 0; i < 3; i++)
683  {
684  Array<OneD, NekDouble> tmp(npts,0.0);
685  Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions+i*npts, 1);
686  }
687 }
688 
689 /**
690  *
691  */
693  const LibUtilities::SessionReaderSharedPtr& pSession,
695 {
696  m_movingBodyCalls = 0;
697  m_session->LoadParameter("Kinvis",m_kinvis);
698  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
699 
700  LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
701  int colrank = vcomm->GetColumnComm()->GetRank();
702  int nproc = vcomm->GetColumnComm()->GetSize();
703 
704  //number of structral modes
705  int npts;
706  bool homostrip;
707  m_session->MatchSolverInfo("HomoStrip","True",homostrip,false);
708 
709  if(!homostrip) //full resolutions
710  {
711  npts = m_session->GetParameter("HomModesZ");
712  }
713  else
714  {
715  m_session->LoadParameter("HomStructModesZ", npts);
716  }
717 
720  m_MotionVars[1] = Array<OneD, NekDouble>(3*npts,0.0);
721 
723  ZIDs = pFields[0]->GetZIDs();
724  m_np = ZIDs.num_elements();
725 
726  std::string vibtype = m_session->GetSolverInfo("VibrationType");
727 
728  if(boost::iequals(vibtype, "Constrained"))
729  {
730  m_vdim = 1;
731  }
732  else if (boost::iequals(vibtype, "Free"))
733  {
734  m_vdim = 2;
735  }
736  else if (boost::iequals(vibtype, "Forced"))
737  {
738  m_vdim = 0;
739  return;
740  }
741 
742  if(!homostrip)
743  {
744  m_session->LoadParameter("LZ", m_lhom);
745  int nplanes = m_session->GetParameter("HomModesZ");
746  m_FFT =
748  "NekFFTW", nplanes);
749  }
750  else
751  {
752  int nstrips;
753  NekDouble DistStrip;
754 
755  m_session->LoadParameter("DistStrip", DistStrip);
756  m_session->LoadParameter("Strip_Z", nstrips);
757  m_lhom = nstrips * DistStrip;
758  m_FFT =
760  "NekFFTW", nstrips);
761  }
762 
763  // load the structural dynamic parameters from xml file
764  m_session->LoadParameter("StructRho", m_structrho);
765  m_session->LoadParameter("StructDamp", m_structdamp, 0.0);
766 
767  // Identify whether the fictitious mass method is active for explicit
768  // coupling of fluid solver and structural dynamics solver
769  bool fictmass;
770  m_session->MatchSolverInfo("FictitiousMassMethod", "True",
771  fictmass, false);
772  if(fictmass)
773  {
774  NekDouble fictrho, fictdamp;
775  m_session->LoadParameter("FictMass", fictrho);
776  m_session->LoadParameter("FictDamp", fictdamp);
777  m_structrho += fictrho;
778  m_structdamp += fictdamp;
779 
780  // Storage array of Struct Velocity and Acceleration used for
781  // extrapolation of fictitious force
784  for (int i = 0; i < m_motion.num_elements(); ++i)
785  {
787  m_fA[i] = Array<OneD, Array<OneD, NekDouble> > (2);
788 
789  for(int n = 0; n < 2; ++n)
790  {
791  m_fV[i][n] = Array<OneD, NekDouble>(npts, 0.0);
792  m_fA[i][n] = Array<OneD, NekDouble>(npts, 0.0);
793  }
794  }
795  }
796 
797  // Setting the coefficient matrices for solving structural dynamic ODEs
798  SetDynEqCoeffMatrix(pFields);
799 
800  // Set initial condition for cable's motion
801  int cnt = 0;
802 
803  for(int j = 0; j < m_funcName.num_elements(); j++)
804  {
805  // loading from the specified files through inputstream
806  if(m_IsFromFile[cnt] && m_IsFromFile[cnt+1])
807  {
808  std::ifstream inputStream;
809  int nzpoints;
810 
811  if(homostrip)
812  {
813  m_session->LoadParameter("HomStructModesZ", nzpoints);
814  }
815  else
816  {
817  nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
818  }
819 
820  if (vcomm->GetRank() == 0)
821  {
822  std::string filename = m_session->GetFunctionFilename(
823  m_funcName[j], m_motion[0]);
824 
825  // Open intputstream for cable motions
826  inputStream.open(filename.c_str());
827 
828  // Import the head string from the file
830  for(int n = 0; n < tmp.num_elements(); n++)
831  {
832  inputStream >> tmp[n];
833  }
834 
835  NekDouble time, z_cds;
836  // Import the motion variables from the file
837  for (int n = 0; n < nzpoints; n++)
838  {
839  inputStream >> setprecision(6) >> time;
840  inputStream >> setprecision(6) >> z_cds;
841  inputStream >> setprecision(8) >> m_MotionVars[0][n];
842  inputStream >> setprecision(8) >> m_MotionVars[0][n+nzpoints];
843  inputStream >> setprecision(8) >> m_MotionVars[0][n+2*nzpoints];
844  inputStream >> setprecision(8) >> m_MotionVars[1][n];
845  inputStream >> setprecision(8) >> m_MotionVars[1][n+nzpoints];
846  inputStream >> setprecision(8) >> m_MotionVars[1][n+2*nzpoints];
847  }
848  // Close inputstream for cable motions
849  inputStream.close();
850  }
851  cnt = cnt + 2;
852  }
853  else //Evaluate from the functions specified in xml file
854  {
855  if(!homostrip)
856  {
857  int nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
858  Array<OneD, NekDouble> z_coords(nzpoints,0.0);
860  = pFields[0]->GetHomogeneousBasis()->GetZ();
861 
862  Vmath::Smul(nzpoints,m_lhom/2.0,pts,1,z_coords,1);
863  Vmath::Sadd(nzpoints,m_lhom/2.0,z_coords,1,z_coords,1);
864 
865  Array<OneD, NekDouble> x0(m_np, 0.0);
866  Array<OneD, NekDouble> x1(m_np, 0.0);
867  Array<OneD, NekDouble> x2(m_np, 0.0);
868  for (int plane = 0; plane < m_np; plane++)
869  {
870  x2[plane] = z_coords[ZIDs[plane]];
871  }
872 
873  Array<OneD, NekDouble> tmp (m_np,0.0);
874  Array<OneD, NekDouble> tmp0(m_np,0.0);
875  Array<OneD, NekDouble> tmp1(m_np,0.0);
876  LibUtilities::EquationSharedPtr ffunc0,ffunc1;
877 
878  ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
879  ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
880  ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
881  ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
882 
883  int offset = j*npts;
884  Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0]+offset,1);
885  Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1]+offset,1);
886 
887  if(colrank == 0)
888  {
889  for (int i = 1; i < nproc; ++i)
890  {
891  vcomm->GetColumnComm()->Recv(i, tmp0);
892  vcomm->GetColumnComm()->Recv(i, tmp1);
893  Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0]+offset+i*m_np,1);
894  Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1]+offset+i*m_np,1);
895  }
896  }
897  else
898  {
899  vcomm->GetColumnComm()->Send(0, tmp0);
900  vcomm->GetColumnComm()->Send(0, tmp1);
901  }
902  }
903  else
904  {
905  if(colrank == 0)
906  {
907  int nstrips;
908  m_session->LoadParameter("Strip_Z", nstrips);
909 
910  ASSERTL0(m_session->DefinesSolverInfo("USEFFT"),
911  "Fourier transformation of cable motion is currently "
912  "implemented only for FFTW module.");
913 
914  NekDouble DistStrip;
915  m_session->LoadParameter("DistStrip", DistStrip);
916 
917  Array<OneD, NekDouble> x0(npts, 0.0);
918  Array<OneD, NekDouble> x1(npts, 0.0);
919  Array<OneD, NekDouble> x2(npts, 0.0);
920  Array<OneD, NekDouble> tmp (npts,0.0);
921  Array<OneD, NekDouble> tmp0(npts,0.0);
922  Array<OneD, NekDouble> tmp1(npts,0.0);
923  for (int plane = 0; plane < npts; plane++)
924  {
925  x2[plane] = plane*DistStrip;
926  }
927  LibUtilities::EquationSharedPtr ffunc0,ffunc1;
928  ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
929  ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
930  ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
931  ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
932 
933  int offset = j*npts;
934  Vmath::Vcopy(npts, tmp0, 1, tmp = m_MotionVars[0]+offset,1);
935  Vmath::Vcopy(npts, tmp1, 1, tmp = m_MotionVars[1]+offset,1);
936  }
937  }
938 
939  cnt = cnt + 2;
940  }
941  }
942 }
943 
944 
945 /**
946  *
947  */
950 {
951  int nplanes;
952 
953  bool homostrip;
954  m_session->MatchSolverInfo("HomoStrip","True",homostrip,false);
955 
956  if(!homostrip)
957  {
958  nplanes = m_session->GetParameter("HomModesZ");
959  }
960  else
961  {
962  m_session->LoadParameter("Strip_Z", nplanes);
963  }
964 
967 
968  NekDouble tmp1, tmp2, tmp3;
969  NekDouble tmp4, tmp5, tmp6, tmp7;
970 
971  // load the structural dynamic parameters from xml file
972  NekDouble cabletension;
973  NekDouble bendingstiff;
974  NekDouble structstiff;
975  m_session->LoadParameter("StructStiff", structstiff, 0.0);
976  m_session->LoadParameter("CableTension", cabletension, 0.0);
977  m_session->LoadParameter("BendingStiff", bendingstiff, 0.0);
978 
979  tmp1 = m_timestep * m_timestep;
980  tmp2 = structstiff / m_structrho;
981  tmp3 = m_structdamp / m_structrho;
982  tmp4 = cabletension / m_structrho;
983  tmp5 = bendingstiff / m_structrho;
984 
985  // solve the ODE in the wave space for cable motion to obtain disp, vel and
986  // accel
987 
988  std::string supptype = m_session->GetSolverInfo("SupportType");
989 
990  for(int plane = 0; plane < nplanes; plane++)
991  {
992  int nel = 3;
993  m_CoeffMat_A[plane]
995  m_CoeffMat_B[plane]
997 
998  // Initialised to avoid compiler warnings.
999  unsigned int K = 0;
1000  NekDouble beta = 0.0;
1001 
1002  if (boost::iequals(supptype, "Free-Free"))
1003  {
1004  K = plane/2;
1005  beta = 2.0 * M_PI/m_lhom;
1006  }
1007  else if(boost::iequals(supptype, "Pinned-Pinned"))
1008  {
1009  K = plane+1;
1010  beta = M_PI/m_lhom;
1011  }
1012  else
1013  {
1014  ASSERTL0(false,
1015  "Unrecognized support type for cable's motion");
1016  }
1017 
1018  tmp6 = beta * K;
1019  tmp6 = tmp6 * tmp6;
1020  tmp7 = tmp6 * tmp6;
1021  tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1022 
1023  (*m_CoeffMat_A[plane])(0,0) = tmp7;
1024  (*m_CoeffMat_A[plane])(0,1) = tmp3;
1025  (*m_CoeffMat_A[plane])(0,2) = 1.0;
1026  (*m_CoeffMat_A[plane])(1,0) = 0.0;
1027  (*m_CoeffMat_A[plane])(1,1) = 1.0;
1028  (*m_CoeffMat_A[plane])(1,2) =-m_timestep/2.0;
1029  (*m_CoeffMat_A[plane])(2,0) = 1.0;
1030  (*m_CoeffMat_A[plane])(2,1) = 0.0;
1031  (*m_CoeffMat_A[plane])(2,2) =-tmp1/4.0;
1032  (*m_CoeffMat_B[plane])(0,0) = 0.0;
1033  (*m_CoeffMat_B[plane])(0,1) = 0.0;
1034  (*m_CoeffMat_B[plane])(0,2) = 0.0;
1035  (*m_CoeffMat_B[plane])(1,0) = 0.0;
1036  (*m_CoeffMat_B[plane])(1,1) = 1.0;
1037  (*m_CoeffMat_B[plane])(1,2) = m_timestep/2.0;
1038  (*m_CoeffMat_B[plane])(2,0) = 1.0;
1039  (*m_CoeffMat_B[plane])(2,1) = m_timestep;
1040  (*m_CoeffMat_B[plane])(2,2) = tmp1/4.0;
1041 
1042  m_CoeffMat_A[plane]->Invert();
1043  (*m_CoeffMat_B[plane]) =
1044  (*m_CoeffMat_A[plane]) * (*m_CoeffMat_B[plane]);
1045  }
1046 }
1047 
1048 
1049 /**
1050  * Function to roll time-level storages to the next step layout.
1051  * The stored data associated with the oldest time-level
1052  * (not required anymore) are moved to the top, where they will
1053  * be overwritten as the solution process progresses.
1054  */
1056 {
1057  int nlevels = input.num_elements();
1058 
1060  tmp = input[nlevels-1];
1061 
1062  for(int n = nlevels-1; n > 0; --n)
1063  {
1064  input[n] = input[n-1];
1065  }
1066 
1067  input[0] = tmp;
1068 }
1069 
1070 /**
1071  *
1072  */
1073 void ForcingMovingBody::CheckIsFromFile(const TiXmlElement* pForce)
1074 {
1075 
1078  m_motion[0] = "x";
1079  m_motion[1] = "y";
1080 
1082  // Loading the x-dispalcement (m_zta) and the y-displacement (m_eta)
1083  // Those two variables are bith functions of z and t and the may come
1084  // from an equation (forced vibration) or from another solver which, given
1085  // the aerodynamic forces at the previous step, calculates the
1086  // displacements.
1087 
1088  //Get the body displacement: m_zta and m_eta
1089  const TiXmlElement* funcNameElmt_D
1090  = pForce->FirstChildElement("DISPLACEMENTS");
1091  ASSERTL0(funcNameElmt_D,
1092  "MOVINGBODYFORCE tag has to specify a function name which "
1093  "prescribes the body displacement as d(z,t).");
1094 
1095  m_funcName[0] = funcNameElmt_D->GetText();
1096  ASSERTL0(m_session->DefinesFunction(m_funcName[0]),
1097  "Function '" + m_funcName[0] + "' not defined.");
1098 
1099  //Get the body velocity of movement: d(m_zta)/dt and d(m_eta)/dt
1100  const TiXmlElement* funcNameElmt_V
1101  = pForce->FirstChildElement("VELOCITIES");
1102  ASSERTL0(funcNameElmt_D,
1103  "MOVINGBODYFORCE tag has to specify a function name which "
1104  "prescribes the body velocity of movement as v(z,t).");
1105 
1106  m_funcName[1] = funcNameElmt_V->GetText();
1107  ASSERTL0(m_session->DefinesFunction(m_funcName[1]),
1108  "Function '" + m_funcName[1] + "' not defined.");
1109 
1110 
1111  //Get the body acceleration: dd(m_zta)/ddt and dd(m_eta)/ddt
1112  const TiXmlElement* funcNameElmt_A
1113  = pForce->FirstChildElement("ACCELERATIONS");
1114  ASSERTL0(funcNameElmt_A,
1115  "MOVINGBODYFORCE tag has to specify a function name which "
1116  "prescribes the body acceleration as a(z,t).");
1117 
1118  m_funcName[2] = funcNameElmt_A->GetText();
1119  ASSERTL0(m_session->DefinesFunction(m_funcName[2]),
1120  "Function '" + m_funcName[2] + "' not defined.");
1121 
1123 
1124  // Check Displacement x
1125  vType = m_session->GetFunctionType(m_funcName[0],m_motion[0]);
1127  {
1128  m_IsFromFile[0] = false;
1129  }
1130  else if (vType == LibUtilities::eFunctionTypeFile)
1131  {
1132  m_IsFromFile[0] = true;
1133  }
1134  else
1135  {
1136  ASSERTL0(false, "The displacements in x must be specified via an "
1137  "equation <E> or a file <F>");
1138  }
1139 
1140  // Check Displacement y
1141  vType = m_session->GetFunctionType(m_funcName[0],m_motion[1]);
1143  {
1144  m_IsFromFile[1] = false;
1145  }
1146  else if (vType == LibUtilities::eFunctionTypeFile)
1147  {
1148  m_IsFromFile[1] = true;
1149  }
1150  else
1151  {
1152  ASSERTL0(false, "The displacements in y must be specified via an "
1153  "equation <E> or a file <F>");
1154  }
1155 
1156  // Check Velocity x
1157  vType = m_session->GetFunctionType(m_funcName[1],m_motion[0]);
1159  {
1160  m_IsFromFile[2] = false;
1161  }
1162  else if (vType == LibUtilities::eFunctionTypeFile)
1163  {
1164  m_IsFromFile[2] = true;
1165  }
1166  else
1167  {
1168  ASSERTL0(false, "The velocities in x must be specified via an equation "
1169  "<E> or a file <F>");
1170  }
1171 
1172  // Check Velocity y
1173  vType = m_session->GetFunctionType(m_funcName[1],m_motion[1]);
1175  {
1176  m_IsFromFile[3] = false;
1177  }
1178  else if (vType == LibUtilities::eFunctionTypeFile)
1179  {
1180  m_IsFromFile[3] = true;
1181  }
1182  else
1183  {
1184  ASSERTL0(false, "The velocities in y must be specified via an equation "
1185  "<E> or a file <F>");
1186  }
1187 
1188  // Check Acceleration x
1189  vType = m_session->GetFunctionType(m_funcName[2],m_motion[0]);
1191  {
1192  m_IsFromFile[4] = false;
1193  }
1194  else if (vType == LibUtilities::eFunctionTypeFile)
1195  {
1196  m_IsFromFile[4] = true;
1197  }
1198  else
1199  {
1200  ASSERTL0(false, "The accelerations in x must be specified via an "
1201  "equation <E> or a file <F>");
1202  }
1203 
1204  // Check Acceleration y
1205  vType = m_session->GetFunctionType(m_funcName[2],m_motion[1]);
1207  {
1208  m_IsFromFile[5] = false;
1209  }
1210  else if (vType == LibUtilities::eFunctionTypeFile)
1211  {
1212  m_IsFromFile[5] = true;
1213  }
1214  else
1215  {
1216  ASSERTL0(false, "The accelerations in y must be specified via an "
1217  "equation <E> or a file <F>");
1218  }
1219 }
1220 
1221 /**
1222  *
1223  */
1225  const LibUtilities::SessionReaderSharedPtr& pSession,
1227  const TiXmlElement* pForce)
1228 {
1229  // Get the outputfile name, output frequency and
1230  // the boundary's ID for the cable's wall
1231  std::string typeStr = pForce->Attribute("TYPE");
1232  std::map<std::string, std::string> vParams;
1233 
1234  const TiXmlElement *param = pForce->FirstChildElement("PARAM");
1235  while (param)
1236  {
1237  ASSERTL0(param->Attribute("NAME"),
1238  "Missing attribute 'NAME' for parameter in filter "
1239  + typeStr + "'.");
1240  std::string nameStr = param->Attribute("NAME");
1241 
1242  ASSERTL0(param->GetText(), "Empty value string for param.");
1243  std::string valueStr = param->GetText();
1244 
1245  vParams[nameStr] = valueStr;
1246 
1247  param = param->NextSiblingElement("PARAM");
1248  }
1249 
1250  // Creat the filter for MovingBody, where we performed the calculation of
1251  // fluid forces and write both the aerodynamic forces and motion variables
1252  // into the output files
1254  AllocateSharedPtr(pSession, vParams);
1255 
1256  // Initialise the object of MovingBody filter
1257  m_MovBodyfilter->Initialise(pFields, 0.0);
1258 
1259 }
1260 
1261 }
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
NekDouble m_structrho
mass of the cable per unit length
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
void CheckIsFromFile(const TiXmlElement *pForce)
NekDouble m_structdamp
damping ratio of the cable
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
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > m_MotionVars
storage for the cable's motion(x,y) variables
void SetDynEqCoeffMatrix(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
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
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
STL namespace.
virtual void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
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:151
GlobalMapping::MappingSharedPtr m_mapping
void Newmark_betaSolver(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int m_movingBodyCalls
number of times the movbody have been called
NekDouble m_kinvis
fluid viscous
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
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
int m_vdim
vibration dimension
LibUtilities::NektarFFTSharedPtr m_FFT
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)
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:95
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
NekDouble m_timestep
time step
boost::shared_ptr< Equation > EquationSharedPtr
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
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
FilterMovingBodySharedPtr m_MovBodyfilter
NekDouble m_lhom
length ratio of the cable
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
int m_np
number of planes per processors
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:70
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:266
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