56 ModuleKey ProcessQualityMetric::className =
59 ProcessQualityMetric::create,
60 "add quality metric to field.");
76 if (
m_f->m_comm->TreatAsRankZero())
78 cout <<
"ProcessQualityMetric: Adding quality metric to field"
86 for (
int i = 0; i <
m_f->m_exp[0]->GetExpSize(); ++i)
90 int offset =
m_f->m_exp[0]->GetPhys_Offset(i);
94 ASSERTL0(q.num_elements() == Elmt->GetTotPoints(),
95 "number of points mismatch");
99 m_f->m_exp[0]->FwdTrans_IterPerExp(phys, coeffs);
101 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
102 m_f->m_exp[0]->GetFieldDefinitions();
103 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
105 for (
int i = 0; i < FieldDef.size(); ++i)
107 FieldDef[i]->m_fields.push_back(
"QualityMetric");
108 m_f->m_exp[0]->AppendFieldData(FieldDef[i], FieldData[i]);
111 m_f->m_fielddef = FieldDef;
112 m_f->m_data = FieldData;
122 vector<Array<OneD, NekDouble> > xy;
123 for (
int i = 0; i < geom->GetNumVerts(); i++)
135 for (
int j = 0; j < b[1]->GetNumPoints(); j++)
137 for (
int i = 0; i < b[0]->GetNumPoints(); i++)
139 NekDouble a1 = 0.5 * (1.0 - u[i]), a2 = 0.5 * (1.0 + u[i]);
140 NekDouble b1 = 0.5 * (1.0 - v[j]), b2 = 0.5 * (1.0 + v[j]);
143 dxdz(0, 0) = 0.5 * (-b1 * xy[0][0] + b1 * xy[1][0] +
144 b2 * xy[2][0] - b2 * xy[3][0]);
145 dxdz(1, 0) = 0.5 * (-b1 * xy[0][1] + b1 * xy[1][1] +
146 b2 * xy[2][1] - b2 * xy[3][1]);
148 dxdz(0, 1) = 0.5 * (-a1 * xy[0][0] - a2 * xy[1][0] +
149 a2 * xy[2][0] + a1 * xy[3][0]);
150 dxdz(1, 1) = 0.5 * (-a1 * xy[0][1] - a2 * xy[1][1] +
151 a2 * xy[2][1] + a1 * xy[3][1]);
160 vector<Array<OneD, NekDouble> > xy;
161 for (
int i = 0; i < geom->GetNumVerts(); i++)
173 for (
int i = 0; i < b[0]->GetNumPoints(); i++)
175 for (
int j = 0; j < b[1]->GetNumPoints(); j++)
178 dxdz(0, 0) = -xy[0][0] / 2.0 + xy[1][0] / 2.0;
180 dxdz(0, 1) = -xy[0][0] / 2.0 + xy[2][0] / 2.0;
182 dxdz(1, 0) = -xy[0][1] / 2.0 + xy[1][1] / 2.0;
184 dxdz(1, 1) = -xy[0][1] / 2.0 + xy[2][1] / 2.0;
193 vector<Array<OneD, NekDouble> > xyz;
194 for (
int i = 0; i < geom->GetNumVerts(); i++)
207 for (
int i = 0; i < b[0]->GetNumPoints(); i++)
209 for (
int j = 0; j < b[1]->GetNumPoints(); j++)
211 for (
int k = 0; k < b[2]->GetNumPoints(); k++)
214 dxdz(0, 0) = -xyz[0][0] / 2.0 + xyz[1][0] / 2.0;
216 dxdz(0, 1) = -xyz[0][0] / 2.0 + xyz[2][0] / 2.0;
218 dxdz(0, 2) = -xyz[0][0] / 2.0 + xyz[3][0] / 2.0;
220 dxdz(1, 0) = -xyz[0][1] / 2.0 + xyz[1][1] / 2.0;
222 dxdz(1, 1) = -xyz[0][1] / 2.0 + xyz[2][1] / 2.0;
224 dxdz(1, 2) = -xyz[0][1] / 2.0 + xyz[3][1] / 2.0;
226 dxdz(2, 0) = -xyz[0][2] / 2.0 + xyz[1][2] / 2.0;
228 dxdz(2, 1) = -xyz[0][2] / 2.0 + xyz[2][2] / 2.0;
230 dxdz(2, 2) = -xyz[0][2] / 2.0 + xyz[3][2] / 2.0;
240 vector<Array<OneD, NekDouble> > xyz;
241 for (
int i = 0; i < geom->GetNumVerts(); i++)
254 for (
int k = 0; k < b[2]->GetNumPoints(); k++)
256 for (
int j = 0; j < b[1]->GetNumPoints(); j++)
258 for (
int i = 0; i < b[0]->GetNumPoints(); i++)
260 NekDouble xi1 = 0.5 * (1 + eta1[i]) * (1 - eta3[k]) - 1.0;
263 b2 = 0.5 * (1 + eta2[j]);
265 c2 = 0.5 * (1 + eta3[k]);
269 dxdz(0, 0) = 0.5 * (-b1 * xyz[0][0] + b1 * xyz[1][0] +
270 b2 * xyz[2][0] - b2 * xyz[3][0]);
271 dxdz(1, 0) = 0.5 * (-b1 * xyz[0][1] + b1 * xyz[1][1] +
272 b2 * xyz[2][1] - b2 * xyz[3][1]);
273 dxdz(2, 0) = 0.5 * (-b1 * xyz[0][2] + b1 * xyz[1][2] +
274 b2 * xyz[2][2] - b2 * xyz[3][2]);
276 dxdz(0, 1) = 0.5 * ((a2 - c1) * xyz[0][0] - a2 * xyz[1][0] +
277 a2 * xyz[2][0] + (c1 - a2) * xyz[3][0] -
278 c2 * xyz[4][0] + c2 * xyz[5][0]);
279 dxdz(1, 1) = 0.5 * ((a2 - c1) * xyz[0][1] - a2 * xyz[1][1] +
280 a2 * xyz[2][1] + (c1 - a2) * xyz[3][1] -
281 c2 * xyz[4][1] + c2 * xyz[5][1]);
282 dxdz(2, 1) = 0.5 * ((a2 - c1) * xyz[0][2] - a2 * xyz[1][2] +
283 a2 * xyz[2][2] + (c1 - a2) * xyz[3][2] -
284 c2 * xyz[4][2] + c2 * xyz[5][2]);
286 dxdz(0, 2) = 0.5 * (-b1 * xyz[0][0] - b2 * xyz[3][0] +
287 b1 * xyz[4][0] + b2 * xyz[5][0]);
288 dxdz(1, 2) = 0.5 * (-b1 * xyz[0][1] - b2 * xyz[3][1] +
289 b1 * xyz[4][1] + b2 * xyz[5][1]);
290 dxdz(2, 2) = 0.5 * (-b1 * xyz[0][2] - b2 * xyz[3][2] +
291 b1 * xyz[4][2] + b2 * xyz[5][2]);
301 vector<Array<OneD, NekDouble> > xyz;
302 for (
int i = 0; i < geom->GetNumVerts(); i++)
315 for (
int k = 0; k < b[2]->GetNumPoints(); k++)
317 for (
int j = 0; j < b[1]->GetNumPoints(); j++)
319 for (
int i = 0; i < b[0]->GetNumPoints(); i++)
324 b2 = 0.5 * (1 + eta2[j]);
326 c2 = 0.5 * (1 + eta3[k]);
331 -0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0] +
332 0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0] -
333 0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0] +
334 0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
336 -0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1] +
337 0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1] -
338 0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1] +
339 0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
341 -0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2] +
342 0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2] -
343 0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2] +
344 0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
347 -0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0] +
348 0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0] -
349 0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0] +
350 0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
352 -0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1] +
353 0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1] -
354 0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1] +
355 0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
357 -0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2] +
358 0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2] -
359 0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2] +
360 0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
363 -0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0] -
364 0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0] +
365 0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0] +
366 0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
368 -0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1] -
369 0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1] +
370 0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1] +
371 0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
373 -0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2] -
374 0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2] +
375 0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2] +
376 0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
401 const int expDim = chi->GetNumBases();
404 vector<LibUtilities::BasisKey> basisKeys;
405 bool needsInterp =
false;
407 for (
int i = 0; i < expDim; ++i)
409 nElemPts *= pElem[i].GetNumPoints();
411 needsInterp || pElem[i].GetNumPoints() < p[i].GetNumPoints() - 1;
417 err <<
"Interpolating from higher order geometry to lower order in "
418 <<
"element " << geom->GetGlobalID();
422 for (
int i = 0; i < expDim; ++i)
426 ? chi->GetBasis(i)->GetBasisKey()
428 chi->GetBasisNumModes(i), pElem[i]));
432 switch (chi->DetShapeType())
436 basisKeys[0], basisKeys[1]);
440 basisKeys[0], basisKeys[1]);
444 basisKeys[0], basisKeys[1], basisKeys[2]);
448 basisKeys[0], basisKeys[1], basisKeys[2]);
452 basisKeys[0], basisKeys[1], basisKeys[2]);
460 const int pts = deriv[0][0].num_elements();
461 const int nq = chiMod->GetTotPoints();
468 for (
int k = 0; k < pts; ++k)
473 for (
int i = 0; i < expDim; ++i)
475 for (
int j = 0; j < expDim; ++j)
477 jac(j, i) = deriv[i][j][k];
481 jacIdeal = jac * i2rm[k];
486 jacDet = jacIdeal(0, 0) * jacIdeal(1, 1) -
487 jacIdeal(0, 1) * jacIdeal(1, 0);
489 else if (expDim == 3)
491 jacDet = jacIdeal(0, 0) * (jacIdeal(1, 1) * jacIdeal(2, 2) -
492 jacIdeal(2, 1) * jacIdeal(1, 2)) -
493 jacIdeal(0, 1) * (jacIdeal(1, 0) * jacIdeal(2, 2) -
494 jacIdeal(2, 0) * jacIdeal(1, 2)) +
495 jacIdeal(0, 2) * (jacIdeal(1, 0) * jacIdeal(2, 1) -
496 jacIdeal(2, 0) * jacIdeal(1, 1));
511 for (
int i = 0; i < expDim; ++i)
513 for (
int j = 0; j < expDim; ++j)
515 frob += jacIdeal(i,j) * jacIdeal(i,j);
519 NekDouble sigma = 0.5*(jacDet + sqrt(jacDet*jacDet));
520 eta[k] = expDim * pow(sigma, 2.0/expDim) / frob;
526 NekDouble mx = -1.0 * numeric_limits<double>::max();
527 NekDouble mn = numeric_limits<double>::max();
528 for(
int k = 0; k < pts; k++)
533 for(
int k = 0; k < pts; k++)
540 if (needsInterp && pts != 1)
548 else if (expDim == 3)
555 ASSERTL0(
false,
"mesh dim makes no sense");
virtual ~ProcessQualityMetric()
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
vector< DNekMat > MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom, StdRegions::StdExpansionSharedPtr chi)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
std::vector< PointsKey > PointsKeyVector
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Represents a command-line configuration option.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
virtual void Process(po::variables_map &vm)
Write mesh to output file.
pair< ModuleType, string > ModuleKey
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
boost::shared_ptr< Field > FieldSharedPtr
boost::shared_ptr< Expansion > ExpansionSharedPtr
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis...
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Abstract base class for processing modules.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Array< OneD, NekDouble > GetQ(LocalRegions::ExpansionSharedPtr e, bool s)
boost::shared_ptr< Geometry > GeometrySharedPtr
Describes the specification for a Basis.
boost::shared_ptr< PointGeom > PointGeomSharedPtr
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.