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