Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FilterMovingBody.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterMovingBody.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: Output moving body forces during time-stepping.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <iomanip>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 
48 std::string FilterMovingBody::className = SolverUtils::GetFilterFactory().
49  RegisterCreatorFunction("MovingBody",
50  FilterMovingBody::create,
51  "Moving Body Filter");
52 /**
53  *
54  */
55 FilterMovingBody::FilterMovingBody(
57  const ParamMap &pParams)
58  : Filter(pSession),
59  m_session(pSession)
60 {
61  ParamMap::const_iterator it;
62 
63  // OutputFile
64  it = pParams.find("OutputFile");
65  if (it == pParams.end())
66  {
67  m_outputFile_fce = pSession->GetSessionName();
68  m_outputFile_mot = pSession->GetSessionName();
69  }
70  else
71  {
72  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
73 
74  m_outputFile_fce = it->second;
75  m_outputFile_mot = it->second;
76  }
77  if (!(m_outputFile_fce.length() >= 4 &&
78  m_outputFile_fce.substr(m_outputFile_fce.length() - 4) == ".fce"))
79  {
80  m_outputFile_fce += ".fce";
81  }
82  if (!(m_outputFile_mot.length() >= 4 &&
83  m_outputFile_mot.substr(m_outputFile_mot.length() - 4) == ".mot"))
84  {
85  m_outputFile_mot += ".mot";
86  }
87 
88  // OutputFrequency
89  it = pParams.find("OutputFrequency");
90  if (it == pParams.end())
91  {
93  }
94  else
95  {
96  LibUtilities::Equation equ(m_session, it->second);
97  m_outputFrequency = floor(equ.Evaluate());
98  }
99 
100  pSession->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
101  ASSERTL0(m_isHomogeneous1D, "Moving Body implemented just for 3D "
102  "Homogeneous 1D discetisations.");
103 
104  // Boundary (to calculate the forces)
105  it = pParams.find("Boundary");
106  ASSERTL0(it != pParams.end(), "Missing parameter 'Boundary'.");
107  ASSERTL0(it->second.length() > 0, "Missing parameter 'Boundary'.");
108  m_BoundaryString = it->second;
109 }
110 
111 
112 /**
113  *
114  */
116 {
117 
118 }
119 
120 
121 /**
122  *
123  */
126  const NekDouble &time)
127 {
128  m_index_f = 0;
129  m_index_m = 0;
131  // Parse the boundary regions into a list.
132  std::string::size_type FirstInd = m_BoundaryString.find_first_of('[') + 1;
133  std::string::size_type LastInd = m_BoundaryString.find_last_of(']') - 1;
134 
135  ASSERTL0(FirstInd <= LastInd,
136  (std::string("Error reading boundary region definition:") +
137  m_BoundaryString).c_str());
138 
139  std::string IndString = m_BoundaryString.substr(FirstInd,
140  LastInd - FirstInd + 1);
141 
142  bool parseGood = ParseUtils::GenerateSeqVector(IndString.c_str(),
144 
145  ASSERTL0(parseGood && !m_boundaryRegionsIdList.empty(),
146  (std::string("Unable to read boundary regions index "
147  "range for FilterAeroForces: ") + IndString).c_str());
148 
149  // determine what boundary regions need to be considered
150  int cnt;
151 
152  unsigned int numBoundaryRegions
153  = pFields[0]->GetBndConditions().num_elements();
154 
156  numBoundaryRegions, 0);
157 
158  SpatialDomains::BoundaryConditions bcs(m_session,pFields[0]->GetGraph());
159 
161  = bcs.GetBoundaryRegions();
162 
163  SpatialDomains::BoundaryRegionCollection::const_iterator it;
164 
165  for (cnt = 0, it = bregions.begin(); it != bregions.end();
166  ++it, cnt++)
167  {
168  if ( std::find(m_boundaryRegionsIdList.begin(),
169  m_boundaryRegionsIdList.end(), it->first) !=
171  {
172  m_boundaryRegionIsInList[cnt] = 1;
173  }
174  }
175 
176  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
177 
178  if (vComm->GetRank() == 0)
179  {
180  // Open output stream for cable forces
181  m_outputStream[0].open(m_outputFile_fce.c_str());
182  m_outputStream[0] << "#";
183  m_outputStream[0].width(7);
184  m_outputStream[0] << "Time";
185  m_outputStream[0].width(15);
186  m_outputStream[0] << "z";
187  m_outputStream[0].width(15);
188  m_outputStream[0] << "Fx (press)";
189  m_outputStream[0].width(15);
190  m_outputStream[0] << "Fx (visc)";
191  m_outputStream[0].width(15);
192  m_outputStream[0] << "Fx (tot)";
193  m_outputStream[0].width(15);
194  m_outputStream[0] << "Fy (press)";
195  m_outputStream[0].width(15);
196  m_outputStream[0] << "Fy (visc)";
197  m_outputStream[0].width(15);
198  m_outputStream[0] << "Fy (tot)";
199  m_outputStream[0] << endl;
200 
201  // Open output stream for cable motions
202  m_outputStream[1].open(m_outputFile_mot.c_str());
203  m_outputStream[1] << "#";
204  m_outputStream[1].width(7);
205  m_outputStream[1] << "Time";
206  m_outputStream[1].width(15);
207  m_outputStream[1] << "z";
208  m_outputStream[1].width(15);
209  m_outputStream[1] << "Disp_x";
210  m_outputStream[1].width(15);
211  m_outputStream[1] << "Vel_x";
212  m_outputStream[1].width(15);
213  m_outputStream[1] << "Acel_x";
214  m_outputStream[1].width(15);
215  m_outputStream[1] << "Disp_y";
216  m_outputStream[1].width(15);
217  m_outputStream[1] << "Vel_y";
218  m_outputStream[1].width(15);
219  m_outputStream[1] << "Acel_y";
220  m_outputStream[1] << endl;
221  }
222 }
223 
224 
225 /**
226  *
227  */
229  const LibUtilities::SessionReaderSharedPtr &pSession,
231  Array<OneD, NekDouble> &Aeroforces,
232  const NekDouble &time)
233 {
234  int n, cnt, elmtid, nq, offset, boundary;
235  int nt = pFields[0]->GetNpoints();
236  int dim = pFields.num_elements()-1;
237 
239  Array<OneD, int> BoundarytoElmtID;
240  Array<OneD, int> BoundarytoTraceID;
242 
247 
251 
255 
256  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
257  LibUtilities::CommSharedPtr vRowComm = vComm->GetRowComm();
258  LibUtilities::CommSharedPtr vColComm = vComm->GetColumnComm();
259  LibUtilities::CommSharedPtr vCColComm = vColComm->GetColumnComm();
260 
261  // set up storage space for forces on all the planes (z-positions)
262  // on each processors some of them will remain empty as we may
263  // have just few planes per processor
264  int Num_z_pos = pFields[0]->GetHomogeneousBasis()->GetNumModes();
265  Array<OneD, NekDouble> Fx(Num_z_pos,0.0);
266  Array<OneD, NekDouble> Fxp(Num_z_pos,0.0);
267  Array<OneD, NekDouble> Fxv(Num_z_pos,0.0);
268  Array<OneD, NekDouble> Fy(Num_z_pos,0.0);
269  Array<OneD, NekDouble> Fyp(Num_z_pos,0.0);
270  Array<OneD, NekDouble> Fyv(Num_z_pos,0.0);
271 
272 
273  NekDouble rho = (pSession->DefinesParameter("rho"))
274  ? (pSession->GetParameter("rho"))
275  : 1;
276  NekDouble mu = rho*pSession->GetParameter("Kinvis");
277 
278  for(int i = 0; i < pFields.num_elements(); ++i)
279  {
280  pFields[i]->SetWaveSpace(false);
281  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
282  pFields[i]->UpdatePhys());
283  pFields[i]->SetPhysState(true);
284  }
285 
286  // Get the number of local planes on the process and their IDs
287  // to properly locate the forces in the Fx, Fy etc. vectors.
289  ZIDs = pFields[0]->GetZIDs();
290  int local_planes = ZIDs.num_elements();
291 
292  // Homogeneous 1D case Compute forces on all WALL boundaries
293  // This only has to be done on the zero (mean) Fourier mode.
294  for(int plane = 0 ; plane < local_planes; plane++)
295  {
296  pFields[0]->GetPlane(plane)->GetBoundaryToElmtMap(BoundarytoElmtID,
297  BoundarytoTraceID);
298  BndExp = pFields[0]->GetPlane(plane)->GetBndCondExpansions();
300 
301  // loop over the types of boundary conditions
302  for(cnt = n = 0; n < BndExp.num_elements(); ++n)
303  {
304  if(m_boundaryRegionIsInList[n] == 1)
305  {
306  for(int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
307  {
308  // find element of this expansion.
309  elmtid = BoundarytoElmtID[cnt];
310  elmt = pFields[0]->GetPlane(plane)->GetExp(elmtid);
311  nq = elmt->GetTotPoints();
312  offset = pFields[0]->GetPlane(plane)
313  ->GetPhys_Offset(elmtid);
314 
315  // Initialise local arrays for the velocity
316  // gradients size of total number of quadrature
317  // points for each element (hence local).
318  for(int j = 0; j < dim; ++j)
319  {
320  gradU[j] = Array<OneD, NekDouble>(nq,0.0);
321  gradV[j] = Array<OneD, NekDouble>(nq,0.0);
322  gradW[j] = Array<OneD, NekDouble>(nq,0.0);
323  }
324 
325  // identify boundary of element
326  boundary = BoundarytoTraceID[cnt];
327 
328  // Extract fields
329  U = pFields[0]->GetPlane(plane)->GetPhys() + offset;
330  V = pFields[1]->GetPlane(plane)->GetPhys() + offset;
331  P = pFields[3]->GetPlane(plane)->GetPhys() + offset;
332 
333  // compute the gradients
334  elmt->PhysDeriv(U,gradU[0],gradU[1]);
335  elmt->PhysDeriv(V,gradV[0],gradV[1]);
336 
337  // Get face 1D expansion from element expansion
338  bc = BndExp[n]->GetExp(i)->as<LocalRegions::Expansion1D>();
339 
340  // number of points on the boundary
341  int nbc = bc->GetTotPoints();
342 
343  // several vectors for computing the forces
344  Array<OneD, NekDouble> Pb(nbc,0.0);
345 
346  for(int j = 0; j < dim; ++j)
347  {
348  fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
349  fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
350  }
351 
352  Array<OneD, NekDouble> drag_t(nbc,0.0);
353  Array<OneD, NekDouble> lift_t(nbc,0.0);
354  Array<OneD, NekDouble> drag_p(nbc,0.0);
355  Array<OneD, NekDouble> lift_p(nbc,0.0);
356  Array<OneD, NekDouble> temp(nbc,0.0);
357  Array<OneD, NekDouble> temp2(nbc,0.0);
358 
359  // identify boundary of element .
360  boundary = BoundarytoTraceID[cnt];
361 
362  // extraction of the pressure and wss on the
363  // boundary of the element
364  elmt->GetEdgePhysVals(boundary,bc,P,Pb);
365 
366  for(int j = 0; j < dim; ++j)
367  {
368  elmt->GetEdgePhysVals(boundary,bc,gradU[j],fgradU[j]);
369  elmt->GetEdgePhysVals(boundary,bc,gradV[j],fgradV[j]);
370  }
371 
372  //normals of the element
373  const Array<OneD, Array<OneD, NekDouble> > &normals
374  = elmt->GetEdgeNormal(boundary);
375 
376  //
377  // Compute viscous tractive forces on wall from
378  //
379  // t_i = - T_ij * n_j (minus sign for force
380  // exerted BY fluid ON wall),
381  //
382  // where
383  //
384  // T_ij = viscous stress tensor (here in Cartesian
385  // coords)
386  // dU_i dU_j
387  // = RHO * KINVIS * ( ---- + ---- ) .
388  // dx_j dx_i
389 
390  //a) DRAG TERMS
391  //-rho*kinvis*(2*du/dx*nx+(du/dy+dv/dx)*ny)
392 
393  Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, drag_t, 1);
394  Vmath::Vmul(nbc, drag_t, 1, normals[1], 1, drag_t, 1);
395 
396  Vmath::Smul(nbc, 2.0, fgradU[0], 1, fgradU[0], 1);
397  Vmath::Vmul(nbc, fgradU[0], 1, normals[0], 1, temp2, 1);
398  Vmath::Smul(nbc, 0.5, fgradU[0], 1, fgradU[0], 1);
399 
400  Vmath::Vadd(nbc, temp2, 1, drag_t, 1, drag_t, 1);
401  Vmath::Smul(nbc, -mu, drag_t, 1, drag_t, 1);
402 
403  //zero temporary storage vector
404  Vmath::Zero(nbc, temp, 0.0);
405  Vmath::Zero(nbc, temp2, 0.0);
406 
407 
408  //b) LIFT TERMS
409  //-rho*kinvis*(2*dv/dy*ny+(du/dy+dv/dx)*nx)
410 
411  Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, lift_t, 1);
412  Vmath::Vmul(nbc, lift_t, 1, normals[0], 1, lift_t, 1);
413 
414  Vmath::Smul(nbc, 2.0, fgradV[1], 1, fgradV[1], 1);
415  Vmath::Vmul(nbc, fgradV[1], 1, normals[1], 1, temp2, 1);
416  Vmath::Smul(nbc, -0.5, fgradV[1], 1, fgradV[1], 1);
417 
418 
419  Vmath::Vadd(nbc, temp2, 1, lift_t, 1, lift_t, 1);
420  Vmath::Smul(nbc, -mu, lift_t, 1, lift_t, 1);
421 
422  // Compute normal tractive forces on all WALL
423  // boundaries
424 
425  Vmath::Vvtvp(nbc, Pb, 1, normals[0], 1, drag_p, 1,
426  drag_p, 1);
427  Vmath::Vvtvp(nbc, Pb, 1, normals[1], 1, lift_p, 1,
428  lift_p, 1);
429 
430  //integration over the boundary
431  Fxv[ZIDs[plane]] += bc->Integral(drag_t);
432  Fyv[ZIDs[plane]] += bc->Integral(lift_t);
433 
434  Fxp[ZIDs[plane]] += bc->Integral(drag_p);
435  Fyp[ZIDs[plane]] += bc->Integral(lift_p);
436  }
437  }
438  else
439  {
440  cnt += BndExp[n]->GetExpSize();
441  }
442  }
443  }
444 
445  for(int i = 0; i < pFields.num_elements(); ++i)
446  {
447  pFields[i]->SetWaveSpace(true);
448  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
449  pFields[i]->UpdatePhys());
450  pFields[i]->SetPhysState(false);
451  }
452 
453  // In case we are using an hybrid parallelisation we need
454  // to reduce the forces on the same plane coming from
455  // different mesh partitions.
456  // It is quite an expensive communication, therefore
457  // we check to make sure it is actually required.
458 
459  if(vComm->GetRowComm()->GetSize() > 0)
460  {
461  // NOTE 1: We can eventually sum the viscous and pressure
462  // component before doing the communication, thus
463  // reducing by a factor 2 the communication.
464  // NOTE 2: We may want to set up in the Comm class an AllReduce
465  // routine wich can handle arrays more efficiently
466  for(int plane = 0 ; plane < local_planes; plane++)
467  {
468  vRowComm->AllReduce(Fxp[ZIDs[plane]], LibUtilities::ReduceSum);
469  vRowComm->AllReduce(Fxv[ZIDs[plane]], LibUtilities::ReduceSum);
470  vRowComm->AllReduce(Fyp[ZIDs[plane]], LibUtilities::ReduceSum);
471  vRowComm->AllReduce(Fyv[ZIDs[plane]], LibUtilities::ReduceSum);
472  }
473  }
474 
475  // At this point on rank (0,n) of the Mesh partion communicator we have
476  // the total areo forces on the planes which are on the same column
477  // communicator. Since the planes are scattered on different processors
478  // some of the entries of the vector Fxp, Fxp etc. are still zero.
479  // Now we need to reduce the values on a single vector on rank (0,0) of the
480  // global communicator.
481  if(!pSession->DefinesSolverInfo("HomoStrip"))
482  {
483  if(vComm->GetRowComm()->GetRank() == 0)
484  {
485  for(int z = 0 ; z < Num_z_pos; z++)
486  {
487  vColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
488  vColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
489  vColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
490  vColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
491  }
492  }
493  }
494  else
495  {
496  if(vComm->GetRowComm()->GetRank() == 0)
497  {
498  for(int z = 0 ; z < Num_z_pos; z++)
499  {
500  vCColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
501  vCColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
502  vCColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
503  vCColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
504  }
505  }
506  }
507 
508  if(!pSession->DefinesSolverInfo("HomoStrip"))
509  {
510  //set the forces imparted on the cable's wall
511  for(int plane = 0 ; plane < local_planes; plane++)
512  {
513  Aeroforces[plane] = Fxp[ZIDs[plane]]
514  + Fxv[ZIDs[plane]];
515  Aeroforces[plane + local_planes] = Fyp[ZIDs[plane]]
516  + Fyv[ZIDs[plane]];
517  }
518 
519  // Only output every m_outputFrequency.
520  if ((m_index_f++) % m_outputFrequency)
521  {
522  return;
523  }
524 
525  // At thi point in rank (0,0) we have the full vectors
526  // containing Fxp,Fxv,Fyp and Fyv where different positions
527  // in the vectors correspond to different planes.
528  // Here we write it to file. We do it just on one porcess
529 
530  Array<OneD, NekDouble> z_coords(Num_z_pos,0.0);
532  = pFields[0]->GetHomogeneousBasis()->GetZ();
533 
534  NekDouble LZ;
535  pSession->LoadParameter("LZ", LZ);
536  Vmath::Smul(Num_z_pos,LZ/2.0,pts,1,z_coords,1);
537  Vmath::Sadd(Num_z_pos,LZ/2.0,z_coords,1,z_coords,1);
538  if (vComm->GetRank() == 0)
539  {
540 
541  Vmath::Vadd(Num_z_pos,Fxp,1,Fxv,1,Fx,1);
542  Vmath::Vadd(Num_z_pos,Fyp,1,Fyv,1,Fy,1);
543 
544  for(int i = 0 ; i < Num_z_pos; i++)
545  {
546  m_outputStream[0].width(8);
547  m_outputStream[0] << setprecision(6) << time;
548 
549  m_outputStream[0].width(15);
550  m_outputStream[0] << setprecision(6) << z_coords[i];
551 
552  m_outputStream[0].width(15);
553  m_outputStream[0] << setprecision(8) << Fxp[i];
554  m_outputStream[0].width(15);
555  m_outputStream[0] << setprecision(8) << Fxv[i];
556  m_outputStream[0].width(15);
557  m_outputStream[0] << setprecision(8) << Fx[i];
558 
559  m_outputStream[0].width(15);
560  m_outputStream[0] << setprecision(8) << Fyp[i];
561  m_outputStream[0].width(15);
562  m_outputStream[0] << setprecision(8) << Fyv[i];
563  m_outputStream[0].width(15);
564  m_outputStream[0] << setprecision(8) << Fy[i];
565  m_outputStream[0] << endl;
566  }
567  }
568  }
569  else
570  {
571  // average the forces over local planes for each processor
572  Array<OneD, NekDouble> fces(6,0.0);
573  for(int plane = 0 ; plane < local_planes; plane++)
574  {
575  fces[0] += Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
576  fces[1] += Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
577  fces[2] += Fxp[ZIDs[plane]] ;
578  fces[3] += Fyp[ZIDs[plane]] ;
579  fces[4] += Fxv[ZIDs[plane]] ;
580  fces[5] += Fyv[ZIDs[plane]] ;
581  }
582 
583  fces[0] = fces[0]/local_planes;
584  fces[1] = fces[1]/local_planes;
585  fces[2] = fces[2]/local_planes;
586  fces[3] = fces[3]/local_planes;
587  fces[4] = fces[4]/local_planes;
588  fces[5] = fces[5]/local_planes;
589 
590  // average the forces over communicators within each strip
591  vCColComm->AllReduce(fces[0], LibUtilities::ReduceSum);
592  vCColComm->AllReduce(fces[1], LibUtilities::ReduceSum);
593  vCColComm->AllReduce(fces[2], LibUtilities::ReduceSum);
594  vCColComm->AllReduce(fces[3], LibUtilities::ReduceSum);
595  vCColComm->AllReduce(fces[4], LibUtilities::ReduceSum);
596  vCColComm->AllReduce(fces[5], LibUtilities::ReduceSum);
597 
598  int npts = vComm->GetColumnComm()->GetColumnComm()->GetSize();
599 
600  fces[0] = fces[0]/npts;
601  fces[1] = fces[1]/npts;
602  fces[2] = fces[2]/npts;
603  fces[3] = fces[3]/npts;
604  fces[4] = fces[4]/npts;
605  fces[5] = fces[5]/npts;
606 
607  for(int plane = 0 ; plane < local_planes; plane++)
608  {
609  Aeroforces[plane] = fces[0];
610  Aeroforces[plane+local_planes] = fces[1];
611  }
612 
613  // Only output every m_outputFrequency.
615  {
616  return;
617  }
618 
619  int colrank = vColComm->GetRank();
620  int nstrips;
621 
622  NekDouble DistStrip;
623 
624  pSession->LoadParameter("Strip_Z", nstrips);
625  pSession->LoadParameter("DistStrip", DistStrip);
626 
627  Array<OneD, NekDouble> z_coords(nstrips);
628  for(int i = 0; i < nstrips; i++)
629  {
630  z_coords[i] = i * DistStrip;
631  }
632 
633  if(colrank == 0)
634  {
635  m_outputStream[0].width(8);
636  m_outputStream[0] << setprecision(6) << time;
637 
638  m_outputStream[0].width(15);
639  m_outputStream[0] << setprecision(6) << z_coords[0];
640 
641  m_outputStream[0].width(15);
642  m_outputStream[0] << setprecision(8) << fces[2];
643  m_outputStream[0].width(15);
644  m_outputStream[0] << setprecision(8) << fces[4];
645  m_outputStream[0].width(15);
646  m_outputStream[0] << setprecision(8) << fces[0];
647 
648  m_outputStream[0].width(15);
649  m_outputStream[0] << setprecision(8) << fces[3];
650  m_outputStream[0].width(15);
651  m_outputStream[0] << setprecision(8) << fces[5];
652  m_outputStream[0].width(15);
653  m_outputStream[0] << setprecision(8) << fces[1];
654  m_outputStream[0] << endl;
655 
656  for(int i = 1; i < nstrips; i++)
657  {
658  vColComm->Recv(i, fces);
659 
660  m_outputStream[0].width(8);
661  m_outputStream[0] << setprecision(6) << time;
662 
663  m_outputStream[0].width(15);
664  m_outputStream[0] << setprecision(6) << z_coords[i];
665 
666  m_outputStream[0].width(15);
667  m_outputStream[0] << setprecision(8) << fces[2];
668  m_outputStream[0].width(15);
669  m_outputStream[0] << setprecision(8) << fces[4];
670  m_outputStream[0].width(15);
671  m_outputStream[0] << setprecision(8) << fces[0];
672 
673  m_outputStream[0].width(15);
674  m_outputStream[0] << setprecision(8) << fces[3];
675  m_outputStream[0].width(15);
676  m_outputStream[0] << setprecision(8) << fces[5];
677  m_outputStream[0].width(15);
678  m_outputStream[0] << setprecision(8) << fces[1];
679  m_outputStream[0] << endl;
680  }
681  }
682  else
683  {
684  for(int i = 1; i < nstrips; i++)
685  {
686  if(colrank == i)
687  {
688  vColComm->Send(0, fces);
689  }
690  }
691  }
692  }
693 }
694 
695 
696 /**
697  *
698  */
700  const LibUtilities::SessionReaderSharedPtr &pSession,
702  Array<OneD, NekDouble> &MotionVars,
703  const NekDouble &time)
704 {
705  // Only output every m_outputFrequency.
706  if ((m_index_m++) % m_outputFrequency)
707  {
708  return;
709  }
710 
711  int npts;
712  // Length of the cable
713  NekDouble Length;
714 
715  if(!pSession->DefinesSolverInfo("HomoStrip"))
716  {
717  pSession->LoadParameter("LZ", Length);
718  npts = m_session->GetParameter("HomModesZ");
719  }
720  else
721  {
722  pSession->LoadParameter("LC", Length);
723  npts = m_session->GetParameter("HomStructModesZ");
724  }
725 
726  NekDouble z_coords;
727  for(int n = 0; n < npts; n++)
728  {
729  z_coords = Length/npts*n;
730  m_outputStream[1].width(8);
731  m_outputStream[1] << setprecision(6) << time;
732  m_outputStream[1].width(15);
733  m_outputStream[1] << setprecision(6) << z_coords;
734  m_outputStream[1].width(15);
735  m_outputStream[1] << setprecision(8) << MotionVars[n];
736  m_outputStream[1].width(15);
737  m_outputStream[1] << setprecision(8) << MotionVars[npts+n];
738  m_outputStream[1].width(15);
739  m_outputStream[1] << setprecision(8) << MotionVars[2*npts+n];
740  m_outputStream[1].width(15);
741  m_outputStream[1] << setprecision(8) << MotionVars[3*npts+n];
742  m_outputStream[1].width(15);
743  m_outputStream[1] << setprecision(8) << MotionVars[4*npts+n];
744  m_outputStream[1].width(15);
745  m_outputStream[1] << setprecision(8) << MotionVars[5*npts+n];
746  m_outputStream[1] << endl;
747  }
748 }
749 
750 
751 /**
752  *
753  */
756  const NekDouble &time)
757 {
758  if (pFields[0]->GetComm()->GetRank() == 0)
759  {
760  m_outputStream[0].close();
761  m_outputStream[1].close();
762  }
763 }
764 
765 
766 /**
767  *
768  */
770 {
771  return true;
772 }
773 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, std::ofstream > m_outputStream
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
STL namespace.
LibUtilities::SessionReaderSharedPtr m_session
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
void UpdateMotion(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &MotionVars, const NekDouble &time)
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
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
NekDouble Evaluate() const
Definition: Equation.h:102
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:206
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
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
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
void UpdateForce(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Aeroforces, const NekDouble &time)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:227
std::vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void 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