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