54 ProcessWSS::create,
"Computes wall shear stress field.");
60 f->m_writeBndFld =
true;
61 f->m_declareExpansionAsContField =
true;
62 m_f->m_fldToBnd =
false;
64 f->m_declareAsNewField =
true;
75 cout <<
"ProcessWSS: Calculating wall shear stress..." << endl;
78 m_f->m_addNormals =
m_config[
"addnormals"].m_beenSet;
81 string bvalues =
m_config[
"bnd"].as<
string>();
83 if(bvalues.compare(
"All") == 0)
86 BndExp =
m_f->m_exp[0]->GetBndCondExpansions();
88 for(
int i = 0; i < BndExp.num_elements(); ++i)
90 m_f->m_bndRegionsToWrite.push_back(i);
96 m_f->m_bndRegionsToWrite),
"Failed to interpret range string");
99 NekDouble kinvis =
m_f->m_session->GetParameter(
"Kinvis");
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)
106 spacedim +=
m_f->m_fielddef[0]->m_numHomogeneousDir;
109 int nfields =
m_f->m_fielddef[0]->m_fields.size();
110 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)");
114 ASSERTL0(
false,
"Error: wss for a 1D problem cannot "
118 int newfields = spacedim + 1;
119 int nshear = spacedim + 1;
120 int nstress = spacedim*spacedim;
121 int ngrad = spacedim*spacedim;
131 for (
int i = 0; i <
m_f->m_exp.size(); ++i)
133 m_f->m_exp[i]->FillBndCondFromField();
136 m_f->m_exp.resize(nfields + newfields);
138 for(i = 0; i < newfields; ++i)
140 m_f->m_exp[nfields + i] =
m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir, var);
145 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_x");
146 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_y");
147 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_mag");
151 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_x");
152 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_y");
153 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_z");
154 m_f->m_fielddef[0]->m_fields.push_back(
"Shear_mag");
158 for(
int b = 0; b <
m_f->m_bndRegionsToWrite.size(); ++b)
160 int bnd =
m_f->m_bndRegionsToWrite[b];
162 for(i = 0; i < newfields; i++)
164 BndExp[i] =
m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd);
166 for(i = 0; i < spacedim; i++)
168 m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
172 int nqb = BndExp[0]->GetTotPoints();
173 int nqe = BndElmtExp[0]->GetTotPoints();
177 for(i = 0; i < ngrad; ++i)
182 for(i = 0; i < nstress; ++i)
188 for(i = 0; i < nstress; ++i)
193 for(i = 0; i < ngrad; ++i)
198 for(i = 0; i < nshear; ++i)
204 for(i = 0; i < spacedim; ++i)
206 velocity[i] = BndElmtExp[i]->GetPhys();
210 for(i = 0; i < spacedim; ++i)
214 BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
219 BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
226 for(i = 0; i < spacedim; ++i)
228 for(j = 0; j < spacedim; ++j)
231 grad[j*spacedim+i], 1,
232 stress[i*spacedim+j], 1);
235 stress[i*spacedim+j], 1);
240 for(j = 0; j < nstress; ++j)
242 m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j],fstress[j]);
247 m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
249 for(i = 0; i < spacedim; ++i)
256 for(i = 0; i < spacedim; ++i)
258 for(j = 0; j < spacedim; ++j)
267 for(i = 0; i < spacedim; ++i)
273 Vmath::Smul(nqb, -1.0, fshear[nshear-1], 1, fshear[nshear-1], 1);
275 for (i = 0; i < spacedim; i++)
280 BndExp[i]->FwdTrans(fshear[i],
281 BndExp[i]->UpdateCoeffs());
286 for(i = 0; i < spacedim; ++i)
292 Vmath::Vsqrt(nqb, fshear[nshear-1], 1, fshear[nshear-1], 1);
293 BndExp[nshear-1]->FwdTrans(fshear[nshear-1],
294 BndExp[nshear-1]->UpdateCoeffs());
#define ASSERTL0(condition, msg)
pair< ModuleType, string > ModuleKey
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
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
map< string, ConfigOption > m_config
List of configuration values.
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.
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< Field > FieldSharedPtr
Represents a command-line configuration option.
void Zero(int n, T *x, const int incx)
Zero vector.
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()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.