Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessQualityMetric.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessQualityMetric.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Compute quality metric of Roca et al.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessQualityMetric.h"
41 
45 #include <StdRegions/StdTriExp.h>
46 #include <StdRegions/StdQuadExp.h>
47 #include <StdRegions/StdTetExp.h>
48 #include <StdRegions/StdPrismExp.h>
49 
50 namespace Nektar
51 {
52 namespace Utilities
53 {
54 
55 ModuleKey ProcessQualityMetric::className =
57  ModuleKey(eProcessModule, "qualitymetric"),
58  ProcessQualityMetric::create,
59  "add quality metric to field.");
60 
61 ProcessQualityMetric::ProcessQualityMetric(FieldSharedPtr f) :
62  ProcessModule(f)
63 {
64 
65 }
66 
68 {
69 }
70 
71 void ProcessQualityMetric::Process(po::variables_map &vm)
72 {
73  if (m_f->m_verbose)
74  {
75  if(m_f->m_comm->GetRank() == 0)
76  {
77  cout << "ProcessQualityMetric: Adding quality metric to field"
78  << endl;
79  }
80  }
81 
82  Array<OneD, NekDouble> &phys = m_f->m_exp[0]->UpdatePhys();
83  Array<OneD, NekDouble> &coeffs = m_f->m_exp[0]->UpdateCoeffs();
84 
85  for(int i =0; i < m_f->m_exp[0]->GetExpSize(); ++i)
86  {
87  // copy Jacobian into field
88  LocalRegions::ExpansionSharedPtr Elmt = m_f->m_exp[0]->GetExp(i);
89  int offset = m_f->m_exp[0]->GetPhys_Offset(i);
90  Array<OneD, NekDouble> q = GetQ(Elmt);
91  Array<OneD, NekDouble> out = phys + offset;
92 
93  ASSERTL0(q.num_elements() == Elmt->GetTotPoints(), "number of points mismatch");
94  Vmath::Vcopy(q.num_elements(), q, 1, out, 1);
95  }
96 
97  m_f->m_exp[0]->FwdTrans_IterPerExp(phys, coeffs);
98 
99  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
100  = m_f->m_exp[0]->GetFieldDefinitions();
101  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
102 
103  for (int i = 0; i < FieldDef.size(); ++i)
104  {
105  FieldDef[i]->m_fields.push_back("QualityMetric");
106  m_f->m_exp[0]->AppendFieldData(FieldDef[i], FieldData[i]);
107  }
108 
109  m_f->m_fielddef = FieldDef;
110  m_f->m_data = FieldData;
111 }
112 
115 {
116  vector<DNekMat> ret;
117 
118  if(geom->GetShapeType() == LibUtilities::eQuadrilateral)
119  {
120  vector<Array<OneD, NekDouble> > xy;
121  for(int i = 0; i < geom->GetNumVerts(); i++)
122  {
123  Array<OneD, NekDouble> loc(2);
124  SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
125  p->GetCoords(loc);
126  xy.push_back(loc);
127  }
128 
130  Array<OneD, NekDouble> u = b[0]->GetZ();
131  Array<OneD, NekDouble> v = b[1]->GetZ();
132 
133  for(int j = 0; j < b[1]->GetNumPoints(); j++)
134  {
135  for(int i = 0; i < b[0]->GetNumPoints(); i++)
136  {
137  NekDouble a1 = 0.5*(1.0-u[i]), a2 = 0.5*(1.0+u[i]);
138  NekDouble b1 = 0.5*(1.0-v[j]), b2 = 0.5*(1.0+v[j]);
139  DNekMat dxdz(2,2,1.0,eFULL);
140 
141  dxdz(0,0) = 0.5*(-b1*xy[0][0] + b1*xy[1][0] + b2*xy[2][0] - b2*xy[3][0]);
142  dxdz(1,0) = 0.5*(-b1*xy[0][1] + b1*xy[1][1] + b2*xy[2][1] - b2*xy[3][1]);
143 
144  dxdz(0,1) = 0.5*(-a1*xy[0][0] - a2*xy[1][0] + a2*xy[2][0] + a1*xy[3][0]);
145  dxdz(1,1) = 0.5*(-a1*xy[0][1] - a2*xy[1][1] + a2*xy[2][1] + a1*xy[3][1]);
146 
147  dxdz.Invert();
148  ret.push_back(dxdz);
149  }
150  }
151  }
152  else if(geom->GetShapeType() == LibUtilities::eTriangle)
153  {
154  vector<Array<OneD, NekDouble> > xy;
155  for(int i = 0; i < geom->GetNumVerts(); i++)
156  {
157  Array<OneD, NekDouble> loc(2);
158  SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
159  p->GetCoords(loc);
160  xy.push_back(loc);
161  }
162 
164  Array<OneD, NekDouble> u = b[0]->GetZ();
165  Array<OneD, NekDouble> v = b[1]->GetZ();
166 
167  for(int i = 0; i < b[0]->GetNumPoints(); i++)
168  {
169  for(int j = 0; j < b[1]->GetNumPoints(); j++)
170  {
171  DNekMat dxdz(2,2,1.0,eFULL);
172  dxdz(0,0) = -xy[0][0]/2.0 + xy[1][0]/2.0;
173 
174  dxdz(0,1) = -xy[0][0]/2.0 + xy[2][0]/2.0;
175 
176  dxdz(1,0) = -xy[0][1]/2.0 + xy[1][1]/2.0;
177 
178  dxdz(1,1) = -xy[0][1]/2.0 + xy[2][1]/2.0;
179 
180  dxdz.Invert();
181  ret.push_back(dxdz);
182  }
183  }
184  }
185  else if(geom->GetShapeType() == LibUtilities::eTetrahedron)
186  {
187  vector<Array<OneD, NekDouble> > xyz;
188  for(int i = 0; i < geom->GetNumVerts(); i++)
189  {
190  Array<OneD, NekDouble> loc(3);
191  SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
192  p->GetCoords(loc);
193  xyz.push_back(loc);
194  }
195 
197  Array<OneD, NekDouble> u = b[0]->GetZ();
198  Array<OneD, NekDouble> v = b[1]->GetZ();
199  Array<OneD, NekDouble> z = b[2]->GetZ();
200 
201  for(int i = 0; i < b[0]->GetNumPoints(); i++)
202  {
203  for(int j = 0; j < b[1]->GetNumPoints(); j++)
204  {
205  for(int k = 0; k < b[2]->GetNumPoints(); k++)
206  {
207  DNekMat dxdz(3,3,1.0,eFULL);
208  dxdz(0,0) = -xyz[0][0]/2.0 + xyz[1][0]/2.0;
209 
210  dxdz(0,1) = -xyz[0][0]/2.0 + xyz[2][0]/2.0;
211 
212  dxdz(0,2) = -xyz[0][0]/2.0 + xyz[3][0]/2.0;
213 
214 
215  dxdz(1,0) = -xyz[0][1]/2.0 + xyz[1][1]/2.0;
216 
217  dxdz(1,1) = -xyz[0][1]/2.0 + xyz[2][1]/2.0;
218 
219  dxdz(1,2) = -xyz[0][1]/2.0 + xyz[3][1]/2.0;
220 
221 
222  dxdz(2,0) = -xyz[0][2]/2.0 + xyz[1][2]/2.0;
223 
224  dxdz(2,1) = -xyz[0][2]/2.0 + xyz[2][2]/2.0;
225 
226  dxdz(2,2) = -xyz[0][2]/2.0 + xyz[3][2]/2.0;
227 
228  dxdz.Invert();
229  ret.push_back(dxdz);
230  }
231  }
232  }
233  }
234  else if(geom->GetShapeType() == LibUtilities::ePrism)
235  {
236  vector<Array<OneD, NekDouble> > xyz;
237  for(int i = 0; i < geom->GetNumVerts(); i++)
238  {
239  Array<OneD, NekDouble> loc(3);
240  SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
241  p->GetCoords(loc);
242  xyz.push_back(loc);
243  }
244 
246  Array<OneD, NekDouble> eta1 = b[0]->GetZ();
247  Array<OneD, NekDouble> eta2 = b[1]->GetZ();
248  Array<OneD, NekDouble> eta3 = b[2]->GetZ();
249 
250  for(int k = 0; k < b[2]->GetNumPoints(); k++)
251  {
252  for(int j = 0; j < b[1]->GetNumPoints(); j++)
253  {
254  for(int i = 0; i < b[0]->GetNumPoints(); i++)
255  {
256  NekDouble xi1 = 0.5*(1+eta1[i])*(1-eta3[k])-1.0;
257  NekDouble a2 = 0.5*(1+xi1);
258  NekDouble b1 = 0.5*(1-eta2[j]), b2 = 0.5*(1+eta2[j]);
259  NekDouble c1 = 0.5*(1-eta3[k]), c2 = 0.5*(1+eta3[k]);
260 
261  DNekMat dxdz(3,3,1.0,eFULL);
262 
263  dxdz(0,0) = 0.5*(-b1*xyz[0][0] + b1*xyz[1][0] + b2*xyz[2][0] - b2*xyz[3][0]);
264  dxdz(1,0) = 0.5*(-b1*xyz[0][1] + b1*xyz[1][1] + b2*xyz[2][1] - b2*xyz[3][1]);
265  dxdz(2,0) = 0.5*(-b1*xyz[0][2] + b1*xyz[1][2] + b2*xyz[2][2] - b2*xyz[3][2]);
266 
267  dxdz(0,1) = 0.5*((a2-c1)*xyz[0][0] - a2*xyz[1][0] + a2*xyz[2][0] + (c1-a2)*xyz[3][0] - c2*xyz[4][0] + c2*xyz[5][0]);
268  dxdz(1,1) = 0.5*((a2-c1)*xyz[0][1] - a2*xyz[1][1] + a2*xyz[2][1] + (c1-a2)*xyz[3][1] - c2*xyz[4][1] + c2*xyz[5][1]);
269  dxdz(2,1) = 0.5*((a2-c1)*xyz[0][2] - a2*xyz[1][2] + a2*xyz[2][2] + (c1-a2)*xyz[3][2] - c2*xyz[4][2] + c2*xyz[5][2]);
270 
271  dxdz(0,2) = 0.5*(-b1*xyz[0][0] - b2*xyz[3][0] + b1*xyz[4][0] + b2*xyz[5][0]);
272  dxdz(1,2) = 0.5*(-b1*xyz[0][1] - b2*xyz[3][1] + b1*xyz[4][1] + b2*xyz[5][1]);
273  dxdz(2,2) = 0.5*(-b1*xyz[0][2] - b2*xyz[3][2] + b1*xyz[4][2] + b2*xyz[5][2]);
274 
275  dxdz.Invert();
276  ret.push_back(dxdz);
277  }
278  }
279  }
280  }
281  else
282  {
283  ASSERTL0(false,"not coded");
284  }
285 
286 
287 
288  return ret;
289 }
290 
292 {
293  SpatialDomains::GeometrySharedPtr geom = e->GetGeom();
294  StdRegions::StdExpansionSharedPtr chi = e->GetGeom()->GetXmap();
295  LibUtilities::PointsKeyVector p = chi->GetPointsKeys();
296  LibUtilities::PointsKeyVector pElem= e->GetPointsKeys();
297  SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
298  const int expDim = chi->GetNumBases();
299  int nElemPts = 1;
300 
301  vector<LibUtilities::BasisKey> basisKeys;
302  bool needsInterp = false;
303 
304  for (int i = 0; i < expDim; ++i)
305  {
306  nElemPts *= pElem[i].GetNumPoints();
307  needsInterp =
308  needsInterp || pElem[i].GetNumPoints() < p[i].GetNumPoints() -1;
309  }
310 
311  if (needsInterp)
312  {
313  stringstream err;
314  err << "Interpolating from higher order geometry to lower order in "
315  << "element " << geom->GetGlobalID();
316  NEKERROR(ErrorUtil::ewarning, err.str());
317  }
318 
319  for (int i = 0; i < expDim; ++i)
320  {
321  basisKeys.push_back(
322  needsInterp ? chi->GetBasis(i)->GetBasisKey() :
323  LibUtilities::BasisKey(chi->GetBasisType(i),
324  chi->GetBasisNumModes(i),
325  pElem[i]));
326  }
327 
329  switch(chi->DetShapeType())
330  {
333  basisKeys[0], basisKeys[1]);
334  break;
337  basisKeys[0], basisKeys[1]);
338  break;
341  basisKeys[0], basisKeys[1], basisKeys[2]);
342  break;
345  basisKeys[0], basisKeys[1], basisKeys[2]);
346  break;
347  default:
348  ASSERTL0(false, "nope");
349  }
350 
351  SpatialDomains::DerivStorage deriv = gfac->GetDeriv(pElem);
352 
353  const int pts = deriv[0][0].num_elements();
354  const int nq = chiMod->GetTotPoints();
355 
356  ASSERTL0(pts == nq, "what");
357 
358  vector<DNekMat> i2rm = MappingIdealToRef(geom, chiMod);
359  Array<OneD, NekDouble> eta(nq);
360 
361  for (int k = 0; k < pts; ++k)
362  {
363  DNekMat jac (expDim, expDim, 0.0, eFULL);
364  DNekMat jacIdeal(expDim, expDim, 0.0, eFULL);
365 
366  for (int i = 0; i < expDim; ++i)
367  {
368  for (int j = 0; j < expDim; ++j)
369  {
370  jac(j,i) = deriv[i][j][k];
371  }
372  }
373 
374  jacIdeal = jac * i2rm[k];
375  NekDouble jacDet;
376 
377  if(expDim == 2)
378  {
379  jacDet = jacIdeal(0,0) * jacIdeal(1,1) - jacIdeal(0,1)*jacIdeal(1,0);
380  }
381  else if(expDim == 3)
382  {
383  jacDet = jacIdeal(0,0) * (jacIdeal(1,1)*jacIdeal(2,2) - jacIdeal(2,1)*jacIdeal(1,2)) -
384  jacIdeal(0,1) * (jacIdeal(1,0)*jacIdeal(2,2) - jacIdeal(2,0)*jacIdeal(1,2)) +
385  jacIdeal(0,2) * (jacIdeal(1,0)*jacIdeal(2,1) - jacIdeal(2,0)*jacIdeal(1,1));
386  }
387  else
388  {
389  ASSERTL0(false,"silly exp dim");
390  }
391 
392  NekDouble frob = 0.0;
393 
394  for (int i = 0; i < expDim; ++i)
395  {
396  for (int j = 0; j < expDim; ++j)
397  {
398  frob += jacIdeal(i,j) * jacIdeal(i,j);
399  }
400  }
401 
402  NekDouble sigma = 0.5*(jacDet + sqrt(jacDet*jacDet));
403  eta[k] = expDim * pow(sigma, 2.0/expDim) / frob;
404  }
405 
406  // Project onto output stuff
407  if (needsInterp && pts != 1)
408  {
409  Array<OneD, NekDouble> tmp(nElemPts);
410 
411  if (expDim == 2)
412  {
413  LibUtilities::Interp2D(p[0], p[1], eta, pElem[0], pElem[1], tmp);
414  }
415  else if(expDim == 3)
416  {
417  LibUtilities::Interp3D(p[0], p[1], p[2], eta, pElem[0], pElem[1],
418  pElem[2], tmp);
419  }
420  else
421  {
422  ASSERTL0(false,"mesh dim makes no sense");
423  }
424 
425  eta = tmp;
426  }
427 
428  if (pts == 1)
429  {
430  Vmath::Fill(nq-1, eta[0], &eta[1], 1);
431  }
432 
433  return eta;
434 }
435 
436 }
437 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:185
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
virtual void Process()=0
STL namespace.
FieldSharedPtr m_f
Field object.
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...
Definition: Interp.cpp:116
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
vector< DNekMat > MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom, StdRegions::StdExpansionSharedPtr chi)
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...
Definition: Interp.cpp:186
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Array< OneD, NekDouble > GetQ(LocalRegions::ExpansionSharedPtr e)
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
Describes the specification for a Basis.
Definition: Basis.h:50
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215