Write mesh to output file.
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)
88 Array<OneD, const MultiRegions::ExpListSharedPtr> BndExp =
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;
129 Array<OneD, Array<OneD, NekDouble> > velocity(nfields), grad(ngrad),
131 Array<OneD, Array<OneD, NekDouble> > stress(nstress), fstress(nstress);
132 Array<OneD, Array<OneD, NekDouble> > fshear(nshear);
134 Array<OneD, MultiRegions::ExpListSharedPtr> BndExp(newfields);
135 Array<OneD, MultiRegions::ExpListSharedPtr> BndElmtExp(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");
166 SpatialDomains::BoundaryConditions bcs(
m_f->m_session,
167 m_f->m_exp[0]->GetGraph());
169 bcs.GetBoundaryRegions();
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)
206 grad[i] = Array<OneD, NekDouble>(nqe);
209 for (i = 0; i < nstress; ++i)
211 stress[i] = Array<OneD, NekDouble>(nqe);
215 for (i = 0; i < nstress; ++i)
217 fstress[i] = Array<OneD, NekDouble>(nqb);
220 for (i = 0; i < ngrad; ++i)
222 fgrad[i] = Array<OneD, NekDouble>(nqb);
225 for (i = 0; i < nshear; ++i)
227 fshear[i] = Array<OneD, NekDouble>(nqb, 0.0);
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]);
274 Array<OneD, Array<OneD, NekDouble> > normals;
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)
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
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
void Neg(int n, T *x, const int incx)
Negate x = -x.
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.
FieldSharedPtr m_f
Field object.