54 "Computes wall shear stress field.");
61 f->m_writeBndFld =
true;
62 f->m_declareExpansionAsContField =
true;
63 f->m_requireBoundaryExpansion =
true;
64 m_f->m_fldToBnd =
false;
75 if (
m_f->m_comm->TreatAsRankZero())
77 cout <<
"ProcessWSS: Calculating wall shear stress..." << endl;
81 m_f->m_addNormals =
m_config[
"addnormals"].m_beenSet;
84 string bvalues =
m_config[
"bnd"].as<
string>();
86 if (bvalues.compare(
"All") == 0)
89 m_f->m_exp[0]->GetBndCondExpansions();
91 for (
int i = 0; i < BndExp.num_elements(); ++i)
93 m_f->m_bndRegionsToWrite.push_back(i);
99 m_f->m_bndRegionsToWrite),
100 "Failed to interpret range string");
103 NekDouble kinvis =
m_f->m_session->GetParameter(
"Kinvis");
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)
110 spacedim +=
m_f->m_fielddef[0]->m_numHomogeneousDir;
113 int nfields =
m_f->m_fielddef[0]->m_fields.size();
115 "Implicit assumption that input is in incompressible format of "
116 "(u,v,p) or (u,v,w,p)");
120 ASSERTL0(
false,
"Error: wss for a 1D problem cannot "
124 int newfields = spacedim + 1;
125 int nshear = spacedim + 1;
126 int nstress = spacedim * spacedim;
127 int ngrad = spacedim * spacedim;
138 for (
int i = 0; i <
m_f->m_exp.size(); ++i)
140 m_f->m_exp[i]->FillBndCondFromField();
143 m_f->m_exp.resize(nfields + newfields);
145 for (i = 0; i < newfields; ++i)
147 m_f->m_exp[nfields + i] =
148 m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir, var);
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");
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");
167 m_f->m_exp[0]->GetGraph());
170 SpatialDomains::BoundaryRegionCollection::const_iterator breg_it;
171 map<int, int> BndRegionMap;
173 for (breg_it = bregions.begin(); breg_it != bregions.end();
176 BndRegionMap[breg_it->first] = cnt;
179 for (
int b = 0; b <
m_f->m_bndRegionsToWrite.size(); ++b)
181 if (BndRegionMap.count(
m_f->m_bndRegionsToWrite[b]) == 1)
183 int bnd = BndRegionMap[
m_f->m_bndRegionsToWrite[b]];
186 for (i = 0; i < newfields; i++)
189 m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd);
191 for (i = 0; i < spacedim; i++)
193 m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
197 int nqb = BndExp[0]->GetTotPoints();
198 int nqe = BndElmtExp[0]->GetTotPoints();
204 for (i = 0; i < ngrad; ++i)
209 for (i = 0; i < nstress; ++i)
215 for (i = 0; i < nstress; ++i)
220 for (i = 0; i < ngrad; ++i)
225 for (i = 0; i < nshear; ++i)
231 for (i = 0; i < spacedim; ++i)
233 velocity[i] = BndElmtExp[i]->GetPhys();
237 for (i = 0; i < spacedim; ++i)
241 BndElmtExp[i]->PhysDeriv(velocity[i],
242 grad[i * spacedim + 0],
243 grad[i * spacedim + 1]);
247 BndElmtExp[i]->PhysDeriv(
248 velocity[i], grad[i * spacedim + 0],
249 grad[i * spacedim + 1], grad[i * spacedim + 2]);
254 for (i = 0; i < spacedim; ++i)
256 for (j = 0; j < spacedim; ++j)
259 grad[j * spacedim + i], 1,
260 stress[i * spacedim + j], 1);
262 Vmath::Smul(nqe, kinvis, stress[i * spacedim + j], 1,
263 stress[i * spacedim + j], 1);
268 for (j = 0; j < nstress; ++j)
270 m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j], fstress[j]);
275 m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
277 for (i = 0; i < spacedim; ++i)
284 for (i = 0; i < spacedim; ++i)
286 for (j = 0; j < spacedim; ++j)
288 Vmath::Vvtvp(nqb, normals[j], 1, fstress[i * spacedim + j],
289 1, fshear[i], 1, fshear[i], 1);
294 for (i = 0; i < spacedim; ++i)
297 fshear[nshear - 1], 1, fshear[nshear - 1], 1);
299 Vmath::Smul(nqb, -1.0, fshear[nshear - 1], 1, fshear[nshear - 1],
302 for (i = 0; i < spacedim; i++)
305 fshear[i], 1, fshear[i], 1);
306 BndExp[i]->FwdTrans(fshear[i], BndExp[i]->UpdateCoeffs());
311 for (i = 0; i < spacedim; ++i)
314 fshear[nshear - 1], 1, fshear[nshear - 1], 1);
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());
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
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)
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
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.
virtual void Process(po::variables_map &vm)
Write mesh to output file.
boost::shared_ptr< Field > FieldSharedPtr
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
void Neg(int n, T *x, const int incx)
Negate x = -x.
const BoundaryRegionCollection & GetBoundaryRegions(void) const
void Zero(int n, T *x, const int incx)
Zero vector.
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.
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.