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