Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
Nektar::Utilities::ProcessWSS Class Reference

This processing module calculates the wall shear stress and adds it as an extra-field to the output file, and writes it to a surface output file. More...

#include <ProcessWSS.h>

Inheritance diagram for Nektar::Utilities::ProcessWSS:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessWSS:
Collaboration graph
[legend]

Public Member Functions

 ProcessWSS (FieldSharedPtr f)
 
virtual ~ProcessWSS ()
 
virtual void Process (po::variables_map &vm)
 Write mesh to output file. More...
 
- Public Member Functions inherited from Nektar::Utilities::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
 ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
 
void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
void PrintConfig ()
 Print out all configuration options for a module. More...
 
void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
bool GetRequireEquiSpaced (void)
 
void SetRequireEquiSpaced (bool pVal)
 
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 
 Module (MeshSharedPtr p_m)
 
virtual void Process ()=0
 
void RegisterConfig (string key, string value)
 
void PrintConfig ()
 
void SetDefaults ()
 
MeshSharedPtr GetMesh ()
 
virtual void ProcessVertices ()
 Extract element vertices. More...
 

Static Public Member Functions

static boost::shared_ptr< Modulecreate (FieldSharedPtr f)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::Utilities::Module
 Module ()
 
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual void ProcessElements ()
 Generate element IDs. More...
 
virtual void ProcessComposites ()
 Generate composites. More...
 
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
void PrismLines (int prism, PerMap &perFaces, set< int > &prismsDone, vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
MeshSharedPtr m_mesh
 Mesh object. More...
 

Detailed Description

This processing module calculates the wall shear stress and adds it as an extra-field to the output file, and writes it to a surface output file.

Definition at line 50 of file ProcessWSS.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessWSS::ProcessWSS ( FieldSharedPtr  f)

Definition at line 56 of file ProcessWSS.cpp.

References Nektar::Utilities::Module::m_config, and Nektar::Utilities::Module::m_f.

56  : 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 }
map< string, ConfigOption > m_config
List of configuration values.
FieldSharedPtr m_f
Field object.
Nektar::Utilities::ProcessWSS::~ProcessWSS ( )
virtual

Definition at line 65 of file ProcessWSS.cpp.

66 {
67 }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessWSS::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 54 of file ProcessWSS.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

54  {
56  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Nektar::Utilities::ProcessWSS::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 69 of file ProcessWSS.cpp.

References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::ParseUtils::GenerateOrderedVector(), Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, Vmath::Smul(), Vmath::Svtsvtp(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vsqrt(), Vmath::Vvtvp(), and Vmath::Vvtvvtp().

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 
98  NekDouble m_kinvis;
99  m_kinvis = m_f->m_session->GetParameter("Kinvis");
100 
101  int i, j;
102  int spacedim = m_f->m_graph->GetSpaceDimension();
103  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
104  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
105  {
106  spacedim = 3;
107  }
108 
109  int nfields = m_f->m_fielddef[0]->m_fields.size();
110  ASSERTL0(nfields == spacedim +1,"Implicit assumption that input is in incompressible format of (u,v,p) or (u,v,w,p)");
111 
112  nfields = nfields - 1;
113  if (spacedim == 1)
114  {
115  ASSERTL0(false, "Error: wss for a 1D problem cannot "
116  "be computed");
117  }
118 
119  int newfields = (spacedim == 2)? 3:4;
120  int nshear = (spacedim == 2)? 3:4;
121  int nstress = (spacedim == 2)? 3:6;
122  int ngrad = nfields*nfields;
123 
124  int n, cnt, elmtid, nq, offset, boundary, nfq;
125  int npoints = m_f->m_exp[0]->GetNpoints();
126  Array<OneD, Array<OneD, NekDouble> > velocity(nfields), grad(ngrad), fgrad(ngrad);
127  Array<OneD, Array<OneD, NekDouble> > stress(nstress), fstress(nstress);
128  Array<OneD, Array<OneD, NekDouble> > outfield(newfields), fshear(nshear);
129 
132  Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
134 
135  m_f->m_exp.resize(newfields);
136  string var = "u";
137  for(i = nfields+1; i < newfields; ++i)
138  {
139  m_f->m_exp[i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir, var);
140  }
141 
142  m_f->m_fielddef[0]->m_fields.resize(newfields);
143  if(spacedim == 2)
144  {
145  m_f->m_fielddef[0]->m_fields[0] = "Shear_x";
146  m_f->m_fielddef[0]->m_fields[1] = "Shear_y";
147  m_f->m_fielddef[0]->m_fields[2] = "Shear_mag";
148  }
149  else
150  {
151  m_f->m_fielddef[0]->m_fields[0] = "Shear_x";
152  m_f->m_fielddef[0]->m_fields[1] = "Shear_y";
153  m_f->m_fielddef[0]->m_fields[2] = "Shear_z";
154  m_f->m_fielddef[0]->m_fields[3] = "Shear_mag";
155  }
156 
157 
158  for (i = 0; i < newfields; ++i)
159  {
160  outfield[i] = Array<OneD, NekDouble>(npoints);
161  }
162  for (i = 0; i < nfields; ++i)
163  {
164  velocity[i] = Array<OneD, NekDouble>(npoints);
165  }
166 
167  m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID, BoundarytoTraceID);
168  //get boundary expansions for each field
169  for(int j = 0; j < newfields; ++j)
170  {
171  BndExp[j] = m_f->m_exp[j]->GetBndCondExpansions();
172  }
173 
174  // loop over the types of boundary conditions
175  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
176  {
177  bool doneBnd = false;
178  // identify if boundary has been defined
179  for(int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
180  {
181  if(n == m_f->m_bndRegionsToWrite[b])
182  {
183  doneBnd = true;
184  for(int i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
185  {
186  // find element and face of this expansion.
187  elmtid = BoundarytoElmtID[cnt];
188  elmt = m_f->m_exp[0]->GetExp(elmtid);
189  nq = elmt->GetTotPoints();
190  offset = m_f->m_exp[0]->GetPhys_Offset(elmtid);
191 
192  // Initialise local arrays for the velocity gradients, and stress components
193  // size of total number of quadrature points for each element (hence local).
194  for(int j = 0; j < ngrad; ++j)
195  {
196  grad[j] = Array<OneD, NekDouble>(nq);
197  }
198 
199  for(int j = 0; j < nstress; ++j)
200  {
201  stress[j] = Array<OneD, NekDouble>(nq);
202  }
203 
204  if(nfields == 2)
205  {
206  ASSERTL0(false, "Error: not implemented in 2D.");
207  }
208  else
209  {
210  // Get face 2D expansion from element expansion
211  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
212  nfq = bc->GetTotPoints();
213 
214  //identify boundary of element looking at.
215  boundary = BoundarytoTraceID[cnt];
216 
217  //Get face normals
218  const SpatialDomains::GeomFactorsSharedPtr m_metricinfo = bc->GetMetricInfo();
219 
221  = elmt->GetFaceNormal(boundary);
222 
223  // initialise arrays
224  for(int j = 0; j < nstress; ++j)
225  {
226  fstress[j] = Array<OneD, NekDouble>(nfq);
227  }
228 
229  for(int j = 0; j < nfields*nfields; ++j)
230  {
231  fgrad[j] = Array<OneD, NekDouble>(nfq);
232  }
233 
234  for(int j = 0; j < nshear; ++j)
235  {
236  fshear[j] = Array<OneD, NekDouble>(nfq);
237  }
238 
239 
240  //Extract Velocities
241  for(int j = 0; j < nfields; ++j)
242  {
243  velocity[j] = m_f->m_exp[j]->GetPhys() + offset;
244  }
245 
246  //Compute gradients (velocity correction scheme method)
247  elmt->PhysDeriv(velocity[0],grad[0],grad[1],grad[2]);
248  elmt->PhysDeriv(velocity[1],grad[3],grad[4],grad[5]);
249  elmt->PhysDeriv(velocity[2],grad[6],grad[7],grad[8]);
250 
251  //Compute stress component terms
252  // t_xx = 2.mu.Ux
253  Vmath::Smul (nq,(2*m_kinvis),grad[0],1,stress[0],1);
254  // tyy = 2.mu.Vy
255  Vmath::Smul (nq,(2*m_kinvis),grad[4],1,stress[1],1);
256  // tzz = 2.mu.Wz
257  Vmath::Smul (nq,(2*m_kinvis),grad[8],1,stress[2],1);
258  // txy = mu.(Uy+Vx)
259  Vmath::Vadd (nq,grad[1],1,grad[3],1,stress[3],1);
260  Vmath::Smul (nq,m_kinvis,stress[3],1,stress[3],1);
261  // txz = mu.(Uz+Wx)
262  Vmath::Vadd (nq,grad[2],1,grad[6],1,stress[4],1);
263  Vmath::Smul (nq,m_kinvis,stress[4],1,stress[4],1);
264  // tyz = mu.(Vz+Wy)
265  Vmath::Vadd (nq,grad[5],1,grad[7],1,stress[5],1);
266  Vmath::Smul (nq,m_kinvis,stress[5],1,stress[5],1);
267 
268 
269  // Get face stress values.
270  for(j = 0; j < nstress; ++j)
271  {
272  elmt->GetFacePhysVals(boundary,bc,stress[j],fstress[j]);
273  }
274 
275  //calcuate wss, and update velocity coeffs in the elemental boundary expansion
276  for (j = 0; j< newfields; j++)
277  {
278  outfield[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[j][n]->GetCoeff_Offset(i);
279  }
280 
281  //surface curved
282  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
283  {
284  // Sx
285  Vmath::Vvtvvtp(nfq,normals[0],1,fstress[0],1,
286  normals[1],1,fstress[3],1,fshear[0],1);
287  Vmath::Vvtvp (nfq,normals[2],1,fstress[4],1,fshear[0],1,fshear[0],1);
288 
289  // Sy
290  Vmath::Vvtvvtp(nfq,normals[0],1,fstress[3],1,
291  normals[1],1,fstress[1],1,fshear[1],1);
292  Vmath::Vvtvp (nfq,normals[2],1,fstress[5],1,fshear[1],1,fshear[1],1);
293 
294  // Sz
295  Vmath::Vvtvvtp(nfq,normals[0],1,fstress[4],1,
296  normals[1],1,fstress[5],1,fshear[2],1);
297  Vmath::Vvtvp (nfq,normals[2],1,fstress[2],1,fshear[2],1,fshear[2],1);
298  }
299  else
300  {
301  // Sx
302  Vmath::Svtsvtp(nfq,normals[0][0],fstress[0],1,
303  normals[1][0],fstress[3],1,fshear[0],1);
304  Vmath::Svtvp(nfq,normals[2][0],fstress[4],1,fshear[0],1,fshear[0],1);
305 
306  // Sy
307  Vmath::Svtsvtp(nfq,normals[0][0],fstress[3],1,
308  normals[1][0],fstress[1],1,fshear[1],1);
309  Vmath::Svtvp(nfq,normals[2][0],fstress[5],1,fshear[1],1,fshear[1],1);
310 
311  // Sz
312  Vmath::Svtsvtp(nfq,normals[0][0],fstress[4],1,
313  normals[1][0],fstress[5],1,fshear[2],1);
314  Vmath::Svtvp(nfq,normals[2][0],fstress[2],1,fshear[2],1,fshear[2],1);
315  }
316 
317  // T = T - (T.n)n
318  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
319  {
320  Vmath::Vvtvvtp(nfq,normals[0],1,fshear[0],1,
321  normals[1],1, fshear[1],1,fshear[3],1);
322  Vmath::Vvtvp (nfq,normals[2],1, fshear[2],1,fshear[3],1,fshear[3],1);
323  Vmath::Smul(nfq, -1.0, fshear[3], 1, fshear[3], 1);
324 
325  for (j = 0; j < nfields; j++)
326  {
327  Vmath::Vvtvp(nfq,normals[j], 1, fshear[3], 1, fshear[j], 1, fshear[j], 1);
328  bc->FwdTrans(fshear[j], outfield[j]);
329  }
330  }
331  else
332  {
333  Vmath::Svtsvtp(nfq,normals[0][0],fshear[0],1,
334  normals[1][0],fshear[1],1,fshear[3],1);
335  Vmath::Svtvp(nfq,normals[2][0],fshear[2],1,fshear[3],1,fshear[3],1);
336  Vmath::Smul(nfq, -1.0, fshear[3], 1,fshear[3], 1);
337 
338  for (j = 0; j < nfields; j++)
339  {
340  Vmath::Svtvp(nfq,normals[j][0],fshear[3],1,fshear[j],1,fshear[j],1);
341  bc->FwdTrans(fshear[j], outfield[j]);
342  }
343  }
344 
345  // Tw
346  Vmath::Vvtvvtp(nfq, fshear[0], 1, fshear[0], 1, fshear[1], 1, fshear[1], 1, fshear[3], 1);
347  Vmath::Vvtvp(nfq, fshear[2], 1, fshear[2], 1, fshear[3], 1, fshear[3], 1);
348  Vmath::Vsqrt(nfq, fshear[3], 1, fshear[3], 1);
349  bc->FwdTrans(fshear[3], outfield[3]);
350 
351  }
352  }
353  }
354  }
355  if(doneBnd == false)
356  {
357  cnt += BndExp[0][n]->GetExpSize();
358  }
359  }
360 
361 
362  for(int j = 0; j < newfields; ++j)
363  {
364  for(int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
365  {
366  m_f->m_exp[j]->UpdateBndCondExpansion(m_f->m_bndRegionsToWrite[b]) = BndExp[j][m_f->m_bndRegionsToWrite[b]];
367  }
368  }
369 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:96
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
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
map< string, ConfigOption > m_config
List of configuration values.
FieldSharedPtr m_f
Field object.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
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
double NekDouble
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:523
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:577
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Geometry is curved or has non-constant factors.
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

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessWSS::className
static
Initial value:
=
ProcessWSS::create, "Computes wall shear stress field.")

Definition at line 57 of file ProcessWSS.h.