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