53 ProcessQCriterion::create,
"Computes Q-Criterion.");
68 cout <<
"ProcessQCriterion: Calculating Q Criterion..." << endl;
72 int expdim =
m_f->m_graph->GetMeshDimension();
73 int spacedim = expdim;
74 if ((
m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
75 (
m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
79 int nfields =
m_f->m_fielddef[0]->m_fields.size();
80 if (spacedim == 1 || spacedim == 2)
82 cerr <<
"\n Error: ProcessQCriterion must be computed for a 3D"
83 " (or quasi-3D) case. \n" << endl;
89 int npoints =
m_f->m_exp[0]->GetNpoints();
103 m_f->m_session->LoadParameter(
"Strip_Z",nstrips,1);
105 m_f->m_exp.resize(nfields*nstrips);
107 for (i = 0; i < nfields*nfields; ++i)
112 for (i = 0; i < addfields; ++i)
124 vector<MultiRegions::ExpListSharedPtr> Exp(nstrips*addfields);
126 for(s = 0; s < nstrips; ++s)
128 for (i = 0; i < nfields; ++i)
130 m_f->m_exp[s*nfields+i]->PhysDeriv(
m_f->m_exp[s*nfields+i]->GetPhys(),
138 grad[1 * nfields + 2], 1,
147 grad[2 * nfields + 0], 1,
156 grad[0 * nfields + 1], 1,
173 for (
int k = 0; k < addfields; ++k)
175 Vmath::Smul(npoints, fac, &omega[k][0], 1, &omega[k][0], 1);
183 grad[0 * nfields + 0], 1,
186 grad[1 * nfields + 1], 1,
189 grad[2 * nfields + 2], 1,
201 grad[1 * nfields + 2], 1,
205 &outfield1[0][0], 1);
209 grad[2 * nfields + 0], 1,
213 &outfield2[0][0], 1);
217 grad[0 * nfields + 1], 1,
221 &outfield3[0][0], 1);
225 &outfield2[0][0], 1);
228 &outfield3[0][0], 1);
230 for (
int k = 0; k < addfields; ++k)
233 &outfield3[k][0], 1);
236 Vmath::Vadd(npoints, &outfield3[0][0], 1, &S[0][0], 1, &S[0][0], 1);
237 Vmath::Vsub(npoints, omega[0], 1, S[0], 1, outfield[0], 1);
239 for (
int k = 0; k < addfields; ++k)
246 for (i = 0; i < addfields; ++i)
248 int n = s*addfields + i;
249 Exp[n] =
m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
250 Exp[n]->UpdatePhys() = outfield[i];
251 Exp[n]->FwdTrans(outfield[i],
252 Exp[n]->UpdateCoeffs());
258 for(s = 0; s < nstrips; ++s)
260 for(i = 0; i < addfields; ++i)
262 it =
m_f->m_exp.begin()+s*(nfields+addfields)+nfields+i;
263 m_f->m_exp.insert(it, Exp[s*addfields+i]);
267 vector<string> outname;
268 outname.push_back(
"Q");
270 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
271 =
m_f->m_exp[0]->GetFieldDefinitions();
272 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
274 for(s = 0; s < nstrips; ++s)
276 for (j = 0; j < nfields + addfields; ++j)
278 for (i = 0; i < FieldDef.size()/nstrips; ++i)
280 int n = s * FieldDef.size()/nstrips + i;
284 FieldDef[n]->m_fields.push_back(outname[j-nfields]);
288 FieldDef[n]->m_fields.push_back(
289 m_f->m_fielddef[0]->m_fields[j]);
291 m_f->m_exp[s*(nfields + addfields)+j]->AppendFieldData(FieldDef[n], FieldData[n]);
296 m_f->m_fielddef = FieldDef;
297 m_f->m_data = FieldData;
pair< ModuleType, string > ModuleKey
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.
boost::shared_ptr< Field > FieldSharedPtr
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual ~ProcessQCriterion()
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.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply 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.