Nektar++
FilterAeroForces.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterAeroForces.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 values of aerodynamic forces during time-stepping.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/Memory/NekMemoryManager.hpp>
37 #include <iomanip>
38 #include <LocalRegions/Expansion1D.h>
39 #include <LocalRegions/Expansion2D.h>
40 #include <LocalRegions/Expansion3D.h>
41 #include <SolverUtils/Filters/FilterAeroForces.h>
42 
43 namespace Nektar
44 {
45 namespace SolverUtils
46 {
47 std::string FilterAeroForces::className =
48  GetFilterFactory().RegisterCreatorFunction(
49  "AeroForces", FilterAeroForces::create);
50 
51 /**
52  *
53  */
54 FilterAeroForces::FilterAeroForces(
55  const LibUtilities::SessionReaderSharedPtr &pSession,
56  const ParamMap &pParams) :
57  Filter(pSession)
58 {
59  ParamMap::const_iterator it;
60 
61  // OutputFile
62  it = pParams.find("OutputFile");
63  if (it == pParams.end())
64  {
65  m_outputFile = m_session->GetSessionName();
66  }
67  else
68  {
69  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
70  m_outputFile = it->second;
71  }
72  if (!(m_outputFile.length() >= 4
73  && m_outputFile.substr(m_outputFile.length() - 4) == ".fce"))
74  {
75  m_outputFile += ".fce";
76  }
77 
78  // OutputFrequency
79  it = pParams.find("OutputFrequency");
80  if (it == pParams.end())
81  {
82  m_outputFrequency = 1;
83  }
84  else
85  {
86  LibUtilities::Equation equ(m_session, it->second);
87  m_outputFrequency = floor(equ.Evaluate());
88  }
89 
90  // OutputPlane
91  m_session->MatchSolverInfo("Homogeneous", "1D",
92  m_isHomogeneous1D, false);
93  if(m_isHomogeneous1D)
94  {
95  it = pParams.find("OutputPlane");
96  if (it == pParams.end())
97  {
98  m_outputPlane = 0;
99  }
100  else
101  {
102  LibUtilities::Equation equ(m_session, it->second);
103  m_outputPlane = floor(equ.Evaluate());
104  }
105  }
106 
107  // Boundary (to calculate forces on)
108  it = pParams.find("Boundary");
109  ASSERTL0(it != pParams.end(), "Missing parameter 'Boundary");
110  ASSERTL0(it->second.length() > 0, "Empty parameter 'Boundary'.");
111  m_BoundaryString = it->second;
112 }
113 
114 
115 /**
116  *
117  */
118 FilterAeroForces::~FilterAeroForces()
119 {
120 
121 }
122 
123 
124 /**
125  *
126  */
127 void FilterAeroForces::v_Initialise(
128  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
129  const NekDouble &time)
130 {
131  // Parse the boundary regions into a list.
132  std::string::size_type FirstInd =
133  m_BoundaryString.find_first_of('[') + 1;
134  std::string::size_type LastInd =
135  m_BoundaryString.find_last_of(']') - 1;
136 
137  ASSERTL0(FirstInd <= LastInd,
138  (std::string("Error reading boundary region definition:") +
139  m_BoundaryString).c_str());
140 
141  std::string IndString =
142  m_BoundaryString.substr(FirstInd, LastInd - FirstInd + 1);
143  bool parseGood = ParseUtils::GenerateSeqVector(IndString.c_str(),
144  m_boundaryRegionsIdList);
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  unsigned int numBoundaryRegions =
152  pFields[0]->GetBndConditions().num_elements();
153  m_boundaryRegionIsInList.insert(m_boundaryRegionIsInList.end(),
154  numBoundaryRegions, 0);
155 
156  SpatialDomains::BoundaryConditions bcs(m_session,
157  pFields[0]->GetGraph());
158  const SpatialDomains::BoundaryRegionCollection &bregions =
159  bcs.GetBoundaryRegions();
160  SpatialDomains::BoundaryRegionCollection::const_iterator it;
161 
162  for (cnt = 0, it = bregions.begin(); it != bregions.end();
163  ++it, cnt++)
164  {
165  if ( std::find(m_boundaryRegionsIdList.begin(),
166  m_boundaryRegionsIdList.end(), it->first) !=
167  m_boundaryRegionsIdList.end() )
168  {
169  m_boundaryRegionIsInList[cnt] = 1;
170  }
171  }
172 
173  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
174 
175  if (vComm->GetRank() == 0)
176  {
177  // Open output stream
178  m_outputStream.open(m_outputFile.c_str());
179  m_outputStream << "#";
180  m_outputStream.width(7);
181  m_outputStream << "Time";
182  m_outputStream.width(25);
183  m_outputStream << "Fx (press)";
184  m_outputStream.width(25);
185  m_outputStream << "Fx (visc)";
186  m_outputStream.width(25);
187  m_outputStream << "Fx (tot)";
188  m_outputStream.width(25);
189  m_outputStream << "Fy (press)";
190  m_outputStream.width(25);
191  m_outputStream << "Fy (visc)";
192  m_outputStream.width(25);
193  m_outputStream << "Fy (tot)";
194  m_outputStream.width(25);
195  m_outputStream << "Fz (press)";
196  m_outputStream.width(25);
197  m_outputStream << "Fz (visc)";
198  m_outputStream.width(25);
199  m_outputStream << "Fz (tot)";
200  m_outputStream << endl;
201  }
202 
203  v_Update(pFields, time);
204 }
205 
206 
207 /**
208  *
209  */
210 void FilterAeroForces::v_Update(
211  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
212  const NekDouble &time)
213 {
214  // Only output every m_outputFrequency.
215  if ((m_index++) % m_outputFrequency)
216  {
217  return;
218  }
219 
220  int n, cnt, elmtid, nq, offset, nt, boundary;
221  nt = pFields[0]->GetNpoints();
222  int dim = pFields.num_elements()-1;
223 
224  StdRegions::StdExpansionSharedPtr elmt;
225  Array<OneD, int> BoundarytoElmtID;
226  Array<OneD, int> BoundarytoTraceID;
227  Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
228 
229  Array<OneD, const NekDouble> P(nt);
230  Array<OneD, const NekDouble> U(nt);
231  Array<OneD, const NekDouble> V(nt);
232  Array<OneD, const NekDouble> W(nt);
233 
234  Array<OneD, Array<OneD, NekDouble> > gradU(dim);
235  Array<OneD, Array<OneD, NekDouble> > gradV(dim);
236  Array<OneD, Array<OneD, NekDouble> > gradW(dim);
237 
238  Array<OneD, Array<OneD, NekDouble> > fgradU(dim);
239  Array<OneD, Array<OneD, NekDouble> > fgradV(dim);
240  Array<OneD, Array<OneD, NekDouble> > fgradW(dim);
241 
242  Array<OneD, NekDouble> values;
243  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
244 
245  NekDouble Fx,Fy,Fz,Fxp,Fxv,Fyp,Fyv,Fzp,Fzv;
246 
247  Fxp = 0.0; // x-component of the force due to pressure difference
248  Fxv = 0.0; // x-component of the force due to viscous stress
249  Fx = 0.0; // x-component of the force (total) Fx = Fxp + Fxv (Drag)
250 
251  Fyp = 0.0; // y-component of the force due to pressure difference
252  Fyv = 0.0; // y-component of the force due to viscous stress
253  Fy = 0.0; // y-component of the force (total) Fy = Fyp + Fyv (Lift)
254 
255  Fzp = 0.0; // z-component of the force due to pressure difference
256  Fzv = 0.0; // z-component of the force due to viscous stress
257  Fz = 0.0; // z-component of the force (total) Fz = Fzp + Fzv (Side)
258 
259  NekDouble rho = (m_session->DefinesParameter("rho"))
260  ? (m_session->GetParameter("rho"))
261  : 1;
262  NekDouble mu = rho*m_session->GetParameter("Kinvis");
263 
264  for(int i = 0; i < pFields.num_elements(); ++i)
265  {
266  pFields[i]->SetWaveSpace(false);
267  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
268  pFields[i]->UpdatePhys());
269  pFields[i]->SetPhysState(true);
270  }
271 
272  // Homogeneous 1D case Compute forces on all WALL boundaries
273  // This only has to be done on the zero (mean) Fourier mode.
274  if(m_isHomogeneous1D)
275  {
276  if(vComm->GetColumnComm()->GetRank() == 0)
277  {
278  pFields[0]->GetPlane(0)->GetBoundaryToElmtMap(
279  BoundarytoElmtID,BoundarytoTraceID);
280  BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
281  StdRegions::StdExpansion1DSharedPtr bc;
282 
283  // loop over the types of boundary conditions
284  for(cnt = n = 0; n < BndExp.num_elements(); ++n)
285  {
286  if(m_boundaryRegionIsInList[n] == 1)
287  {
288  for(int i = 0; i < BndExp[n]->GetExpSize();
289  ++i, cnt++)
290  {
291  // find element of this expansion.
292  elmtid = BoundarytoElmtID[cnt];
293  elmt = pFields[0]->GetPlane(0)->GetExp(elmtid);
294  nq = elmt->GetTotPoints();
295  offset = pFields[0]->GetPlane(0)->GetPhys_Offset(elmtid);
296 
297  // Initialise local arrays for the velocity
298  // gradients size of total number of quadrature
299  // points for each element (hence local).
300  for(int j = 0; j < dim; ++j)
301  {
302  gradU[j] = Array<OneD, NekDouble>(nq,0.0);
303  gradV[j] = Array<OneD, NekDouble>(nq,0.0);
304  gradW[j] = Array<OneD, NekDouble>(nq,0.0);
305  }
306 
307  // identify boundary of element
308  boundary = BoundarytoTraceID[cnt];
309 
310  // Extract fields
311  U = pFields[0]->GetPlane(0)->GetPhys() + offset;
312  V = pFields[1]->GetPlane(0)->GetPhys() + offset;
313  P = pFields[3]->GetPlane(0)->GetPhys() + offset;
314 
315  // compute the gradients
316  elmt->PhysDeriv(U,gradU[0],gradU[1]);
317  elmt->PhysDeriv(V,gradV[0],gradV[1]);
318 
319  // Get face 1D expansion from element expansion
320  bc = BndExp[n]->GetExp(i)
321  ->as<LocalRegions::Expansion1D> ();
322 
323  // number of points on the boundary
324  int nbc = bc->GetTotPoints();
325 
326  // several vectors for computing the forces
327  Array<OneD, NekDouble> Pb(nbc,0.0);
328 
329  for(int j = 0; j < dim; ++j)
330  {
331  fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
332  fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
333  }
334 
335  Array<OneD, NekDouble> drag_t(nbc,0.0);
336  Array<OneD, NekDouble> lift_t(nbc,0.0);
337  Array<OneD, NekDouble> drag_p(nbc,0.0);
338  Array<OneD, NekDouble> lift_p(nbc,0.0);
339  Array<OneD, NekDouble> temp(nbc,0.0);
340  Array<OneD, NekDouble> temp2(nbc,0.0);
341 
342  // identify boundary of element .
343  boundary = BoundarytoTraceID[cnt];
344 
345  // extraction of the pressure and wss on the
346  // boundary of the element
347  elmt->GetEdgePhysVals(boundary,bc,P,Pb);
348 
349  for(int j = 0; j < dim; ++j)
350  {
351  elmt->GetEdgePhysVals(boundary,bc,gradU[j],
352  fgradU[j]);
353  elmt->GetEdgePhysVals(boundary,bc,gradV[j],
354  fgradV[j]);
355  }
356 
357  //normals of the element
358  const Array<OneD, Array<OneD, NekDouble> > &normals
359  = elmt->GetEdgeNormal(boundary);
360 
361  //
362  // Compute viscous tractive forces on wall from
363  //
364  // t_i = - T_ij * n_j (minus sign for force
365  // exerted BY fluid ON wall),
366  //
367  // where
368  //
369  // T_ij = viscous stress tensor (here in Cartesian
370  // coords)
371  // dU_i dU_j
372  // = RHO * KINVIS * ( ---- + ---- ) .
373  // dx_j dx_i
374 
375  //a) DRAG TERMS
376  //-rho*kinvis*(2*du/dx*nx+(du/dy+dv/dx)*ny
377 
378  Vmath::Vadd(nbc,fgradU[1],1,fgradV[0],1,drag_t,1);
379  Vmath::Vmul(nbc,drag_t,1,normals[1],1,drag_t,1);
380 
381  Vmath::Smul(nbc,2.0,fgradU[0],1,fgradU[0],1);
382  Vmath::Vmul(nbc,fgradU[0],1,normals[0],1,temp2,1);
383  Vmath::Smul(nbc,0.5,fgradU[0],1,fgradU[0],1);
384 
385  Vmath::Vadd(nbc,temp2,1,drag_t,1,drag_t,1);
386  Vmath::Smul(nbc,-mu,drag_t,1,drag_t,1);
387 
388  //zero temporary storage vector
389  Vmath::Zero(nbc,temp,0);
390  Vmath::Zero(nbc,temp2,0);
391 
392 
393  //b) LIFT TERMS
394  //-rho*kinvis*(2*dv/dy*nx+(du/dy+dv/dx)*nx
395 
396  Vmath::Vadd(nbc,fgradU[1],1,fgradV[0],1,lift_t,1);
397  Vmath::Vmul(nbc,lift_t,1,normals[0],1,lift_t,1);
398 
399  Vmath::Smul(nbc,2.0,fgradV[1],1,fgradV[1],1);
400  Vmath::Vmul(nbc,fgradV[1],1,normals[1],1,temp2,1);
401  Vmath::Smul(nbc,-0.5,fgradV[1],1,fgradV[1],1);
402 
403 
404  Vmath::Vadd(nbc,temp2,1,lift_t,1,lift_t,1);
405  Vmath::Smul(nbc,-mu,lift_t,1,lift_t,1);
406 
407  // Compute normal tractive forces on all WALL
408  // boundaries
409 
410  Vmath::Vvtvp(nbc,Pb,1,normals[0],1,
411  drag_p,1,drag_p, 1);
412  Vmath::Vvtvp(nbc,Pb,1,normals[1],1,
413  lift_p,1,lift_p,1);
414 
415  //integration over the boundary
416  Fxv += bc->Integral(drag_t);
417  Fyv += bc->Integral(lift_t);
418 
419  Fxp += bc->Integral(drag_p);
420  Fyp += bc->Integral(lift_p);
421  }
422  }
423  else
424  {
425  cnt += BndExp[n]->GetExpSize();
426  }
427  }
428  }
429 
430  for(int i = 0; i < pFields.num_elements(); ++i)
431  {
432  pFields[i]->SetWaveSpace(true);
433  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
434  pFields[i]->UpdatePhys());
435  pFields[i]->SetPhysState(false);
436  }
437  }
438  //3D WALL case
439  else if(dim==3 && !m_isHomogeneous1D)
440  {
441  pFields[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
442  BoundarytoTraceID);
443  BndExp = pFields[0]->GetBndCondExpansions();
444  LocalRegions::Expansion2DSharedPtr bc;
445 
446  // loop over the types of boundary conditions
447  for(cnt = n = 0; n < BndExp.num_elements(); ++n)
448  {
449  if(m_boundaryRegionIsInList[n] == 1)
450  {
451  for(int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
452  {
453  // find element of this expansion.
454  elmtid = BoundarytoElmtID[cnt];
455  elmt = pFields[0]->GetExp(elmtid);
456  nq = elmt->GetTotPoints();
457  offset = pFields[0]->GetPhys_Offset(elmtid);
458 
459  // Initialise local arrays for the velocity
460  // gradients size of total number of quadrature
461  // points for each element (hence local).
462  for(int j = 0; j < dim; ++j)
463  {
464  gradU[j] = Array<OneD, NekDouble>(nq,0.0);
465  gradV[j] = Array<OneD, NekDouble>(nq,0.0);
466  gradW[j] = Array<OneD, NekDouble>(nq,0.0);
467  }
468 
469  //identify boundary of element
470  boundary = BoundarytoTraceID[cnt];
471 
472  //Extract fields
473  U = pFields[0]->GetPhys() + offset;
474  V = pFields[1]->GetPhys() + offset;
475  W = pFields[2]->GetPhys() + offset;
476  P = pFields[3]->GetPhys() + offset;
477 
478  //compute the gradients
479  elmt->PhysDeriv(U,gradU[0],gradU[1],gradU[2]);
480  elmt->PhysDeriv(V,gradV[0],gradV[1],gradV[2]);
481  elmt->PhysDeriv(W,gradW[0],gradW[1],gradW[2]);
482 
483  // Get face 2D expansion from element expansion
484  bc = BndExp[n]->GetExp(i)
485  ->as<LocalRegions::Expansion2D> ();
486 
487  //number of points on the boundary
488  int nbc = bc->GetTotPoints();
489 
490  //several vectors for computing the forces
491  Array<OneD, NekDouble> Pb(nbc,0.0);
492 
493  for(int j = 0; j < dim; ++j)
494  {
495  fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
496  fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
497  fgradW[j] = Array<OneD, NekDouble>(nbc,0.0);
498 
499  }
500 
501  Array<OneD, NekDouble> drag_t(nbc,0.0);
502  Array<OneD, NekDouble> lift_t(nbc,0.0);
503  Array<OneD, NekDouble> side_t(nbc,0.0);
504  Array<OneD, NekDouble> drag_p(nbc,0.0);
505  Array<OneD, NekDouble> lift_p(nbc,0.0);
506  Array<OneD, NekDouble> side_p(nbc,0.0);
507  Array<OneD, NekDouble> temp(nbc,0.0);
508  Array<OneD, NekDouble> temp2(nbc,0.0);
509 
510  // identify boundary of element .
511  boundary = BoundarytoTraceID[cnt];
512 
513  // extraction of the pressure and wss on the
514  // boundary of the element
515  elmt->GetFacePhysVals(boundary,bc,P,Pb);
516 
517  for(int j = 0; j < dim; ++j)
518  {
519  elmt->GetFacePhysVals(boundary,bc,gradU[j],
520  fgradU[j]);
521  elmt->GetFacePhysVals(boundary,bc,gradV[j],
522  fgradV[j]);
523  elmt->GetFacePhysVals(boundary,bc,gradW[j],
524  fgradW[j]);
525  }
526 
527  // normals of the element
528  const Array<OneD, Array<OneD, NekDouble> > &normals
529  = elmt->GetFaceNormal(boundary);
530 
531  //
532  // Compute viscous tractive forces on wall from
533  //
534  // t_i = - T_ij * n_j (minus sign for force
535  // exerted BY fluid ON wall),
536  //
537  // where
538  //
539  // T_ij = viscous stress tensor (here in Cartesian
540  // coords)
541  // dU_i dU_j
542  // = RHO * KINVIS * ( ---- + ---- ) .
543  // dx_j dx_i
544 
545  //a) DRAG TERMS
546  //-rho*kinvis*
547  // (2*du/dx*nx+(du/dy+dv/dx)*ny+(du/dz+dw/dx)*nz)
548  Vmath::Vadd(nbc,fgradU[2],1,fgradW[0],1,temp,1);
549  Vmath::Neg(nbc,temp,1);
550  Vmath::Vmul(nbc,temp,1,normals[2],1,temp,1);
551 
552  Vmath::Vadd(nbc,fgradU[1],1,fgradV[0],1,drag_t,1);
553  Vmath::Neg(nbc,drag_t,1);
554  Vmath::Vmul(nbc,drag_t,1,normals[1],1,drag_t,1);
555 
556  Vmath::Smul(nbc,-2.0,fgradU[0],1,fgradU[0],1);
557  Vmath::Vmul(nbc,fgradU[0],1,normals[0],1,temp2,1);
558  Vmath::Smul(nbc,-0.5,fgradU[0],1,fgradU[0],1);
559 
560  Vmath::Vadd(nbc,temp,1,temp2,1,temp,1);
561  Vmath::Vadd(nbc,temp,1,drag_t,1,drag_t,1);
562  Vmath::Smul(nbc,mu,drag_t,1,drag_t,1);
563 
564  //zero temporary storage vector
565  Vmath::Zero(nbc,temp,0);
566  Vmath::Zero(nbc,temp2,0);
567 
568 
569  //b) LIFT TERMS
570  //-rho*kinvis*
571  // (2*dv/dy*nx+(du/dy+dv/dx)*nx+(dv/dz+dw/dy)*nz)
572  Vmath::Vadd(nbc,fgradV[2],1,fgradW[1],1,temp,1);
573  Vmath::Neg(nbc,temp,1);
574  Vmath::Vmul(nbc,temp,1,normals[2],1,temp,1);
575 
576  Vmath::Vadd(nbc,fgradU[1],1,fgradV[0],1,lift_t,1);
577  Vmath::Neg(nbc,lift_t,1);
578  Vmath::Vmul(nbc,lift_t,1,normals[0],1,lift_t,1);
579 
580  Vmath::Smul(nbc,-2.0,fgradV[1],1,fgradV[1],1);
581  Vmath::Vmul(nbc,fgradV[1],1,normals[1],1,temp2,1);
582  Vmath::Smul(nbc,-0.5,fgradV[1],1,fgradV[1],1);
583 
584  Vmath::Vadd(nbc,temp,1,temp2,1,temp,1);
585  Vmath::Vadd(nbc,temp,1,lift_t,1,lift_t,1);
586  Vmath::Smul(nbc,mu,lift_t,1,lift_t,1);
587 
588  //zero temporary storage vector
589  Vmath::Zero(nbc,temp,0);
590  Vmath::Zero(nbc,temp2,0);
591 
592  //b) SIDE TERMS
593  //-rho*kinvis*
594  // (2*dv/dy*nx+(du/dy+dv/dx)*nx+(dv/dz+dw/dy)*nz)
595  Vmath::Vadd(nbc,fgradV[2],1,fgradW[1],1,temp,1);
596  Vmath::Neg(nbc,temp,1);
597  Vmath::Vmul(nbc,temp,1,normals[1],1,temp,1);
598 
599  Vmath::Vadd(nbc,fgradU[2],1,fgradW[0],1,side_t,1);
600  Vmath::Neg(nbc,side_t,1);
601  Vmath::Vmul(nbc,side_t,1,normals[0],1,side_t,1);
602 
603  Vmath::Smul(nbc,-2.0,fgradW[2],1,fgradW[2],1);
604  Vmath::Vmul(nbc,fgradW[2],1,normals[2],1,temp2,1);
605  Vmath::Smul(nbc,-0.5,fgradW[2],1,fgradW[2],1);
606 
607  Vmath::Vadd(nbc,temp,1,temp2,1,temp,1);
608  Vmath::Vadd(nbc,temp,1,side_t,1,side_t,1);
609  Vmath::Smul(nbc,mu,side_t,1,side_t,1);
610 
611 
612  // Compute normal tractive forces on all WALL
613  // boundaries
614  Vmath::Vvtvp(nbc,Pb,1,normals[0],1,
615  drag_p,1,drag_p,1);
616  Vmath::Vvtvp(nbc,Pb,1,normals[1],1,
617  lift_p,1,lift_p,1);
618  Vmath::Vvtvp(nbc,Pb,1,normals[2],1,
619  side_p,1,side_p,1);
620 
621  //integration over the boundary
622  Fxv += bc->Expansion::Integral(drag_t);
623  Fyv += bc->Expansion::Integral(lift_t);
624  Fzv += bc->Expansion::Integral(side_t);
625 
626  Fxp += bc->Expansion::Integral(drag_p);
627  Fyp += bc->Expansion::Integral(lift_p);
628  Fzp += bc->Expansion::Integral(side_p);
629  }
630  }
631  else
632  {
633  cnt += BndExp[n]->GetExpSize();
634  }
635  }
636  }
637  //2D WALL Condition
638  else
639  {
640  pFields[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
641  BoundarytoTraceID);
642  BndExp = pFields[0]->GetBndCondExpansions();
643  StdRegions::StdExpansion1DSharedPtr bc;
644 
645  // loop over the types of boundary conditions
646  for(cnt = n = 0; n < BndExp.num_elements(); ++n)
647  {
648  if(m_boundaryRegionIsInList[n] == 1)
649  {
650  for(int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
651  {
652 
653  elmtid = BoundarytoElmtID[cnt];
654  elmt = pFields[0]->GetExp(elmtid);
655  nq = elmt->GetTotPoints();
656  offset = pFields[0]->GetPhys_Offset(elmtid);
657 
658  for(int j = 0; j < dim; ++j)
659  {
660  gradU[j] = Array<OneD, NekDouble>(nq,0.0);
661  gradV[j] = Array<OneD, NekDouble>(nq,0.0);
662  }
663 
664  boundary = BoundarytoTraceID[cnt];
665 
666  U = pFields[0]->GetPhys() + offset;
667  V = pFields[1]->GetPhys() + offset;
668  P = pFields[2]->GetPhys() + offset;
669 
670  elmt->PhysDeriv(U,gradU[0],gradU[1]);
671  elmt->PhysDeriv(V,gradV[0],gradV[1]);
672 
673  bc = BndExp[n]->GetExp(i)
674  ->as<LocalRegions::Expansion1D> ();
675 
676  int nbc = bc->GetTotPoints();
677  Array<OneD, NekDouble> Pb(nbc,0.0);
678 
679  Array<OneD, NekDouble> drag_t(nbc,0.0);
680  Array<OneD, NekDouble> lift_t(nbc,0.0);
681  Array<OneD, NekDouble> drag_p(nbc,0.0);
682  Array<OneD, NekDouble> lift_p(nbc,0.0);
683  Array<OneD, NekDouble> temp(nbc,0.0);
684 
685  boundary = BoundarytoTraceID[cnt];
686 
687  elmt->GetEdgePhysVals(boundary,bc,P,Pb);
688 
689  for(int j = 0; j < dim; ++j)
690  {
691  fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
692  fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
693 
694  }
695 
696  for(int j = 0; j < dim; ++j)
697  {
698  elmt->GetEdgePhysVals(boundary,bc,gradU[j],
699  fgradU[j]);
700  elmt->GetEdgePhysVals(boundary,bc,gradV[j],
701  fgradV[j]);
702  }
703 
704  const Array<OneD, Array<OneD, NekDouble> > &normals
705  = elmt->GetEdgeNormal(boundary);
706 
707  Vmath::Vadd(nbc,fgradU[1],1,fgradV[0],1,drag_t,1);
708  Vmath::Neg(nbc,drag_t,1);
709  Vmath::Vmul(nbc,drag_t,1,normals[1],1,drag_t,1);
710 
711  Vmath::Smul(nbc,-2.0,fgradU[0],1,fgradU[0],1);
712  Vmath::Vmul(nbc,fgradU[0],1,normals[0],1,temp,1);
713  Vmath::Vadd(nbc,temp,1,drag_t,1,drag_t,1);
714  Vmath::Smul(nbc,mu,drag_t,1,drag_t,1);
715 
716  Vmath::Vadd(nbc,fgradU[1],1,fgradV[0],1,lift_t,1);
717  Vmath::Neg(nbc,lift_t,1);
718  Vmath::Vmul(nbc,lift_t,1,normals[0],1,lift_t,1);
719  Vmath::Smul(nbc,-2.0,fgradV[1],1,fgradV[1],1);
720  Vmath::Vmul(nbc,fgradV[1],1,normals[1],1,temp,1);
721  Vmath::Vadd(nbc,temp,1,lift_t,1,lift_t,1);
722  Vmath::Smul(nbc,mu,lift_t,1,lift_t,1);
723 
724  Vmath::Vvtvp(nbc,Pb,1,normals[0],1,
725  drag_p,1,drag_p,1);
726  Vmath::Vvtvp(nbc,Pb,1,normals[1],1,
727  lift_p,1,lift_p,1);
728 
729  Fxp += bc->Integral(drag_p);
730  Fyp += bc->Integral(lift_p);
731 
732  Fxv += bc->Integral(drag_t);
733  Fyv += bc->Integral(lift_t);
734  }
735  }
736  else
737  {
738  cnt += BndExp[n]->GetExpSize();
739  }
740 
741  }
742 
743  }
744 
745  vComm->AllReduce(Fxp, LibUtilities::ReduceSum);
746  vComm->AllReduce(Fxv, LibUtilities::ReduceSum);
747  Fx = Fxp + Fxv;
748 
749  vComm->AllReduce(Fyp, LibUtilities::ReduceSum);
750  vComm->AllReduce(Fyv, LibUtilities::ReduceSum);
751  Fy = Fyp + Fyv;
752 
753  vComm->AllReduce(Fzp, LibUtilities::ReduceSum);
754  vComm->AllReduce(Fzv, LibUtilities::ReduceSum);
755  Fz = Fzp + Fzv;
756 
757 
758  if (vComm->GetRank() == 0)
759  {
760  m_outputStream.width(8);
761  m_outputStream << setprecision(6) << time;
762 
763  m_outputStream.width(25);
764  m_outputStream << setprecision(8) << Fxp;
765  m_outputStream.width(25);
766  m_outputStream << setprecision(8) << Fxv;
767  m_outputStream.width(25);
768  m_outputStream << setprecision(8) << Fx;
769 
770  m_outputStream.width(25);
771  m_outputStream << setprecision(8) << Fyp;
772  m_outputStream.width(25);
773  m_outputStream << setprecision(8) << Fyv;
774  m_outputStream.width(25);
775  m_outputStream << setprecision(8) << Fy;
776 
777  m_outputStream.width(25);
778  m_outputStream << setprecision(8) << Fzp;
779  m_outputStream.width(25);
780  m_outputStream << setprecision(8) << Fzv;
781  m_outputStream.width(25);
782  m_outputStream << setprecision(8) << Fz;
783 
784  m_outputStream << endl;
785  }
786 }
787 
788 
789 /**
790  *
791  */
792 void FilterAeroForces::v_Finalise(
793  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
794  const NekDouble &time)
795 {
796  if (pFields[0]->GetComm()->GetRank() == 0)
797  {
798  m_outputStream.close();
799  }
800 }
801 
802 
803 /**
804  *
805  */
806 bool FilterAeroForces::v_IsTimeDependent()
807 {
808  return true;
809 }
810 }
811 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
FilterAeroForces.h
Nektar::SolverUtils::FilterAeroForces::FilterAeroForces
SOLVER_UTILS_EXPORT FilterAeroForces(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
Definition: FilterAeroForces.cpp:54
Expansion2D.h
Nektar::SolverUtils::FilterAeroForces::~FilterAeroForces
virtual SOLVER_UTILS_EXPORT ~FilterAeroForces()
Definition: FilterAeroForces.cpp:118
Vmath::Vvtvp
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
Nektar::SolverUtils::FilterAeroForces::v_Update
virtual void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Definition: FilterAeroForces.cpp:210
Nektar::SolverUtils::FilterAeroForces::v_Finalise
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Definition: FilterAeroForces.cpp:792
Nektar::SolverUtils::FilterAeroForces::m_boundaryRegionIsInList
vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
Definition: FilterAeroForces.h:87
Nektar::SolverUtils::FilterAeroForces::m_boundaryRegionsIdList
vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
Definition: FilterAeroForces.h:84
Nektar::SolverUtils::FilterAeroForces::m_outputPlane
unsigned int m_outputPlane
plane to take history point from if using a homogeneous1D expansion
Definition: FilterAeroForces.h:92
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
Nektar::ParseUtils::GenerateSeqVector
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:78
Nektar::SolverUtils::FilterAeroForces::className
static std::string className
Name of the class.
Definition: FilterAeroForces.h:63
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
NekMemoryManager.hpp
Nektar::LibUtilities::Equation::Evaluate
NekDouble Evaluate() const
Definition: Equation.h:102
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:204
Nektar::StdRegions::StdExpansion1DSharedPtr
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
Definition: StdExpansion1D.h:112
Expansion3D.h
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::Filter::ParamMap
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
Nektar::SolverUtils::Filter::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
Nektar::SolverUtils::FilterAeroForces::v_Initialise
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Definition: FilterAeroForces.cpp:127
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:312
Nektar::SolverUtils::GetFilterFactory
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:225
Expansion1D.h
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1791
Nektar::SolverUtils::FilterAeroForces::create
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
Definition: FilterAeroForces.h:52
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
Nektar::LocalRegions::Expansion2DSharedPtr
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
Generated on Tue Dec 8 2015 23:45:28 for Nektar++ by   doxygen 1.8.9.1