Write mesh to output file.
68 if (
m_f->m_comm->TreatAsRankZero())
70 cout <<
"ProcessQCriterion: Calculating Q Criterion..." << endl;
75 int expdim =
m_f->m_graph->GetMeshDimension();
76 int spacedim = expdim;
77 if ((
m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
78 (
m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
82 int nfields =
m_f->m_fielddef[0]->m_fields.size();
83 if (spacedim == 1 || spacedim == 2)
85 cerr <<
"\n Error: ProcessQCriterion must be computed for a 3D"
86 " (or quasi-3D) case. \n"
93 int npoints =
m_f->m_exp[0]->GetNpoints();
95 Array<OneD, Array<OneD, NekDouble> > grad(nfields * nfields);
97 Array<OneD, Array<OneD, NekDouble> > omega(nfields * nfields);
98 Array<OneD, Array<OneD, NekDouble> > S(nfields * nfields);
100 Array<OneD, Array<OneD, NekDouble> > outfield(addfields);
101 Array<OneD, Array<OneD, NekDouble> > outfield1(addfields);
102 Array<OneD, Array<OneD, NekDouble> > outfield2(addfields);
103 Array<OneD, Array<OneD, NekDouble> > outfield3(addfields);
107 m_f->m_session->LoadParameter(
"Strip_Z", nstrips, 1);
109 m_f->m_exp.resize(nfields * nstrips);
111 for (i = 0; i < nfields * nfields; ++i)
113 grad[i] = Array<OneD, NekDouble>(npoints);
116 for (i = 0; i < addfields; ++i)
119 outfield[i] = Array<OneD, NekDouble>(npoints);
120 outfield1[i] = Array<OneD, NekDouble>(npoints);
121 outfield2[i] = Array<OneD, NekDouble>(npoints);
122 outfield3[i] = Array<OneD, NekDouble>(npoints);
124 omega[i] = Array<OneD, NekDouble>(npoints);
125 S[i] = Array<OneD, NekDouble>(npoints);
128 vector<MultiRegions::ExpListSharedPtr> Exp(nstrips * addfields);
130 for (s = 0; s < nstrips; ++s)
132 for (i = 0; i < nfields; ++i)
134 m_f->m_exp[s * nfields + i]->PhysDeriv(
135 m_f->m_exp[s * nfields + i]->GetPhys(), grad[i * nfields],
136 grad[i * nfields + 1], grad[i * nfields + 2]);
140 Vmath::Vsub(npoints, grad[2 * nfields + 1], 1, grad[1 * nfields + 2], 1,
143 Vmath::Vmul(npoints, outfield1[0], 1, outfield1[0], 1, outfield1[0], 1);
146 Vmath::Vsub(npoints, grad[0 * nfields + 2], 1, grad[2 * nfields + 0], 1,
149 Vmath::Vmul(npoints, outfield2[0], 1, outfield2[0], 1, outfield2[0], 1);
152 Vmath::Vsub(npoints, grad[1 * nfields + 0], 1, grad[0 * nfields + 1], 1,
155 Vmath::Vmul(npoints, outfield3[0], 1, outfield3[0], 1, outfield3[0], 1);
160 Vmath::Vadd(npoints, &outfield1[0][0], 1, &outfield2[0][0], 1,
162 Vmath::Vadd(npoints, &omega[0][0], 1, &outfield3[0][0], 1, &omega[0][0],
165 for (
int k = 0; k < addfields; ++k)
167 Vmath::Smul(npoints, fac, &omega[k][0], 1, &omega[k][0], 1);
174 Vmath::Vmul(npoints, grad[0 * nfields + 0], 1, grad[0 * nfields + 0], 1,
176 Vmath::Vmul(npoints, grad[1 * nfields + 1], 1, grad[1 * nfields + 1], 1,
178 Vmath::Vmul(npoints, grad[2 * nfields + 2], 1, grad[2 * nfields + 2], 1,
181 Vmath::Vadd(npoints, &outfield1[0][0], 1, &outfield2[0][0], 1, &S[0][0],
183 Vmath::Vadd(npoints, &S[0][0], 1, &outfield3[0][0], 1, &S[0][0], 1);
186 Vmath::Vadd(npoints, grad[2 * nfields + 1], 1, grad[1 * nfields + 2], 1,
188 Vmath::Vmul(npoints, &outfield1[0][0], 1, &outfield1[0][0], 1,
189 &outfield1[0][0], 1);
192 Vmath::Vadd(npoints, grad[0 * nfields + 2], 1, grad[2 * nfields + 0], 1,
194 Vmath::Vmul(npoints, &outfield2[0][0], 1, &outfield2[0][0], 1,
195 &outfield2[0][0], 1);
198 Vmath::Vadd(npoints, grad[1 * nfields + 0], 1, grad[0 * nfields + 1], 1,
200 Vmath::Vmul(npoints, &outfield3[0][0], 1, &outfield3[0][0], 1,
201 &outfield3[0][0], 1);
203 Vmath::Vadd(npoints, &outfield1[0][0], 1, &outfield2[0][0], 1,
204 &outfield2[0][0], 1);
205 Vmath::Vadd(npoints, &outfield2[0][0], 1, &outfield3[0][0], 1,
206 &outfield3[0][0], 1);
208 for (
int k = 0; k < addfields; ++k)
210 Vmath::Smul(npoints, fac, &outfield3[k][0], 1, &outfield3[k][0], 1);
213 Vmath::Vadd(npoints, &outfield3[0][0], 1, &S[0][0], 1, &S[0][0], 1);
214 Vmath::Vsub(npoints, omega[0], 1, S[0], 1, outfield[0], 1);
216 for (
int k = 0; k < addfields; ++k)
218 Vmath::Smul(npoints, fac, &outfield[k][0], 1, &outfield[k][0], 1);
221 for (i = 0; i < addfields; ++i)
223 int n = s * addfields + i;
225 m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
226 Exp[n]->UpdatePhys() = outfield[i];
227 Exp[n]->FwdTrans(outfield[i], Exp[n]->UpdateCoeffs());
233 for (s = 0; s < nstrips; ++s)
235 for (i = 0; i < addfields; ++i)
237 it =
m_f->m_exp.begin() + s * (nfields + addfields) + nfields + i;
238 m_f->m_exp.insert(it, Exp[s * addfields + i]);
242 vector<string> outname;
243 outname.push_back(
"Q");
245 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
246 m_f->m_exp[0]->GetFieldDefinitions();
247 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
249 for (s = 0; s < nstrips; ++s)
251 for (j = 0; j < nfields + addfields; ++j)
253 for (i = 0; i < FieldDef.size() / nstrips; ++i)
255 int n = s * FieldDef.size() / nstrips + i;
259 FieldDef[n]->m_fields.push_back(outname[j - nfields]);
263 FieldDef[n]->m_fields.push_back(
264 m_f->m_fielddef[0]->m_fields[j]);
266 m_f->m_exp[s * (nfields + addfields) + j]->AppendFieldData(
267 FieldDef[n], FieldData[n]);
272 m_f->m_fielddef = FieldDef;
273 m_f->m_data = FieldData;
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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
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.
FieldSharedPtr m_f
Field object.