54 ProcessWSS::create,
"Computes wall shear stress field.");
60 f->m_writeBndFld =
true;
61 f->m_declareExpansionAsContField =
true;
62 m_f->m_fldToBnd =
false;
73 cout <<
"ProcessWSS: Calculating wall shear stress..." << endl;
76 m_f->m_addNormals =
m_config[
"addnormals"].m_beenSet;
79 string bvalues =
m_config[
"bnd"].as<
string>();
81 if(bvalues.compare(
"All") == 0)
84 BndExp =
m_f->m_exp[0]->GetBndCondExpansions();
86 for(
int i = 0; i < BndExp.num_elements(); ++i)
88 m_f->m_bndRegionsToWrite.push_back(i);
94 m_f->m_bndRegionsToWrite),
"Failed to interpret range string");
97 NekDouble kinvis =
m_f->m_session->GetParameter(
"Kinvis");
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)
104 spacedim +=
m_f->m_fielddef[0]->m_numHomogeneousDir;
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)");
112 ASSERTL0(
false,
"Error: wss for a 1D problem cannot "
116 int newfields = spacedim + 1;
117 int nshear = spacedim + 1;
118 int nstress = spacedim*spacedim;
119 int ngrad = spacedim*spacedim;
129 for (
int i = 0; i <
m_f->m_exp.size(); ++i)
131 m_f->m_exp[i]->FillBndCondFromField();
134 m_f->m_exp.resize(nfields + newfields);
136 for(i = 0; i < newfields; ++i)
138 m_f->m_exp[nfields + i] =
m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir, var);
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");
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");
156 for(
int b = 0; b <
m_f->m_bndRegionsToWrite.size(); ++b)
158 int bnd =
m_f->m_bndRegionsToWrite[b];
160 for(i = 0; i < newfields; i++)
162 BndExp[i] =
m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd);
164 for(i = 0; i < spacedim; i++)
166 m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
170 int nqb = BndExp[0]->GetTotPoints();
171 int nqe = BndElmtExp[0]->GetTotPoints();
175 for(i = 0; i < ngrad; ++i)
180 for(i = 0; i < nstress; ++i)
186 for(i = 0; i < nstress; ++i)
191 for(i = 0; i < ngrad; ++i)
196 for(i = 0; i < nshear; ++i)
202 for(i = 0; i < spacedim; ++i)
204 velocity[i] = BndElmtExp[i]->GetPhys();
208 for(i = 0; i < spacedim; ++i)
212 BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
217 BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
224 for(i = 0; i < spacedim; ++i)
226 for(j = 0; j < spacedim; ++j)
229 grad[j*spacedim+i], 1,
230 stress[i*spacedim+j], 1);
233 stress[i*spacedim+j], 1);
238 for(j = 0; j < nstress; ++j)
240 m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j],fstress[j]);
245 m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
247 for(i = 0; i < spacedim; ++i)
254 for(i = 0; i < spacedim; ++i)
256 for(j = 0; j < spacedim; ++j)
265 for(i = 0; i < spacedim; ++i)
271 Vmath::Smul(nqb, -1.0, fshear[nshear-1], 1, fshear[nshear-1], 1);
273 for (i = 0; i < spacedim; i++)
278 BndExp[i]->FwdTrans(fshear[i],
279 BndExp[i]->UpdateCoeffs());
284 for(i = 0; i < spacedim; ++i)
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());
#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)
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
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.