Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessWSS.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessWSS.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: Computes wss field.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessWSS.h"
41 
44 #include <MultiRegions/ExpList.h>
45 
46 namespace Nektar
47 {
48 namespace Utilities
49 {
50 
51 ModuleKey ProcessWSS::className =
53  ModuleKey(eProcessModule, "wss"),
54  ProcessWSS::create, "Computes wall shear stress field.");
55 
56 ProcessWSS::ProcessWSS(FieldSharedPtr f) : ProcessModule(f)
57 {
58  m_config["bnd"] = ConfigOption(false,"All","Boundary to be extracted");
59  m_config["addnormals"] = ConfigOption(true,"NotSet","Add normals to output");
60  f->m_writeBndFld = true;
61  f->m_declareExpansionAsContField = true;
62  m_f->m_fldToBnd = false;
63 }
64 
66 {
67 }
68 
69 void ProcessWSS::Process(po::variables_map &vm)
70 {
71  if (m_f->m_verbose)
72  {
73  cout << "ProcessWSS: Calculating wall shear stress..." << endl;
74  }
75 
76  m_f->m_addNormals = m_config["addnormals"].m_beenSet;
77 
78  // Set up Field options to output boundary fld
79  string bvalues = m_config["bnd"].as<string>();
80 
81  if(bvalues.compare("All") == 0)
82  {
84  BndExp = m_f->m_exp[0]->GetBndCondExpansions();
85 
86  for(int i = 0; i < BndExp.num_elements(); ++i)
87  {
88  m_f->m_bndRegionsToWrite.push_back(i);
89  }
90  }
91  else
92  {
94  m_f->m_bndRegionsToWrite),"Failed to interpret range string");
95  }
96 
97  NekDouble kinvis = m_f->m_session->GetParameter("Kinvis");
98 
99  int i, j;
100  int spacedim = m_f->m_graph->GetSpaceDimension();
101  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
102  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
103  {
104  spacedim += m_f->m_fielddef[0]->m_numHomogeneousDir;
105  }
106 
107  int nfields = m_f->m_fielddef[0]->m_fields.size();
108  ASSERTL0(m_f->m_fielddef[0]->m_fields[0] == "u","Implicit assumption that input is in incompressible format of (u,v,p) or (u,v,w,p)");
109 
110  if (spacedim == 1)
111  {
112  ASSERTL0(false, "Error: wss for a 1D problem cannot "
113  "be computed");
114  }
115 
116  int newfields = spacedim + 1;
117  int nshear = spacedim + 1;
118  int nstress = spacedim*spacedim;
119  int ngrad = spacedim*spacedim;
120 
121  Array<OneD, Array<OneD, NekDouble> > velocity(nfields), grad(ngrad), fgrad(ngrad);
122  Array<OneD, Array<OneD, NekDouble> > stress(nstress), fstress(nstress);
123  Array<OneD, Array<OneD, NekDouble> > fshear(nshear);
124 
127 
128  // Extract original fields to boundary (for output)
129  for (int i = 0; i < m_f->m_exp.size(); ++i)
130  {
131  m_f->m_exp[i]->FillBndCondFromField();
132  }
133 
134  m_f->m_exp.resize(nfields + newfields);
135  string var = "u";
136  for(i = 0; i < newfields; ++i)
137  {
138  m_f->m_exp[nfields + i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir, var);
139  }
140 
141  if(spacedim == 2)
142  {
143  m_f->m_fielddef[0]->m_fields.push_back("Shear_x");
144  m_f->m_fielddef[0]->m_fields.push_back("Shear_y");
145  m_f->m_fielddef[0]->m_fields.push_back("Shear_mag");
146  }
147  else
148  {
149  m_f->m_fielddef[0]->m_fields.push_back("Shear_x");
150  m_f->m_fielddef[0]->m_fields.push_back("Shear_y");
151  m_f->m_fielddef[0]->m_fields.push_back("Shear_z");
152  m_f->m_fielddef[0]->m_fields.push_back("Shear_mag");
153  }
154 
155  // Loop over boundaries to Write
156  for(int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
157  {
158  int bnd = m_f->m_bndRegionsToWrite[b];
159  // Get expansion list for boundary and for elements containing this bnd
160  for(i = 0; i < newfields; i++)
161  {
162  BndExp[i] = m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd);
163  }
164  for(i = 0; i < spacedim; i++)
165  {
166  m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
167  }
168 
169  // Get number of points in expansions
170  int nqb = BndExp[0]->GetTotPoints();
171  int nqe = BndElmtExp[0]->GetTotPoints();
172 
173  // Initialise local arrays for the velocity gradients, and stress components
174  // size of total number of quadrature points for elements in this bnd
175  for(i = 0; i < ngrad; ++i)
176  {
177  grad[i] = Array<OneD, NekDouble>(nqe);
178  }
179 
180  for(i = 0; i < nstress; ++i)
181  {
182  stress[i] = Array<OneD, NekDouble>(nqe);
183  }
184 
185  // initialise arrays in the boundary
186  for(i = 0; i < nstress; ++i)
187  {
188  fstress[i] = Array<OneD, NekDouble>(nqb);
189  }
190 
191  for(i = 0; i < ngrad; ++i)
192  {
193  fgrad[i] = Array<OneD, NekDouble>(nqb);
194  }
195 
196  for(i = 0; i < nshear; ++i)
197  {
198  fshear[i] = Array<OneD, NekDouble>(nqb, 0.0);
199  }
200 
201  //Extract Velocities
202  for(i = 0; i < spacedim; ++i)
203  {
204  velocity[i] = BndElmtExp[i]->GetPhys();
205  }
206 
207  //Compute gradients (velocity correction scheme method)
208  for(i = 0; i < spacedim; ++i)
209  {
210  if (spacedim == 2)
211  {
212  BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
213  grad[i*spacedim+1]);
214  }
215  else
216  {
217  BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
218  grad[i*spacedim+1],
219  grad[i*spacedim+2]);
220  }
221  }
222 
223  //Compute stress component terms tau_ij = mu*(u_i,j + u_j,i)
224  for(i = 0; i < spacedim; ++i)
225  {
226  for(j = 0; j < spacedim; ++j)
227  {
228  Vmath::Vadd(nqe, grad[i*spacedim+j], 1,
229  grad[j*spacedim+i], 1,
230  stress[i*spacedim+j], 1);
231 
232  Vmath::Smul(nqe, kinvis, stress[i*spacedim+j], 1,
233  stress[i*spacedim+j], 1);
234  }
235  }
236 
237  // Get boundary stress values.
238  for(j = 0; j < nstress; ++j)
239  {
240  m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j],fstress[j]);
241  }
242 
243  //Get normals
245  m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
246  // Reverse normals, to get correct orientation for the body
247  for(i = 0; i < spacedim; ++i)
248  {
249  Vmath::Neg(nqb, normals[i], 1);
250  }
251 
252  //calculate wss, and update coeffs in the boundary expansion
253  // S = tau_ij * n_j
254  for(i = 0; i < spacedim; ++i)
255  {
256  for(j = 0; j < spacedim; ++j)
257  {
258  Vmath::Vvtvp(nqb,normals[j],1,fstress[i*spacedim+j],1,
259  fshear[i],1,
260  fshear[i],1);
261  }
262  }
263 
264  // T = S - (S.n)n
265  for(i = 0; i < spacedim; ++i)
266  {
267  Vmath::Vvtvp(nqb,normals[i],1,fshear[i],1,
268  fshear[nshear-1],1,
269  fshear[nshear-1],1);
270  }
271  Vmath::Smul(nqb, -1.0, fshear[nshear-1], 1, fshear[nshear-1], 1);
272 
273  for (i = 0; i < spacedim; i++)
274  {
275  Vmath::Vvtvp(nqb,normals[i], 1, fshear[nshear-1], 1,
276  fshear[i], 1,
277  fshear[i], 1);
278  BndExp[i]->FwdTrans(fshear[i],
279  BndExp[i]->UpdateCoeffs());
280  }
281 
282  // Tw
283  Vmath::Zero(nqb, fshear[nshear-1], 1);
284  for(i = 0; i < spacedim; ++i)
285  {
286  Vmath::Vvtvp(nqb,fshear[i],1,fshear[i],1,
287  fshear[nshear-1],1,
288  fshear[nshear-1],1);
289  }
290  Vmath::Vsqrt(nqb, fshear[nshear-1], 1, fshear[nshear-1], 1);
291  BndExp[nshear-1]->FwdTrans(fshear[nshear-1],
292  BndExp[nshear-1]->UpdateCoeffs());
293  }
294 
295 }
296 
297 }
298 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:97
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
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.
FieldSharedPtr m_f
Field object.
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
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
Represents a command-line configuration option.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215