Nektar++
Points.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Points.hpp
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Header file of Points definition
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #ifndef NEKTAR_LIB_UTILITIES_FOUNDATIONS_POINTS_H
36 #define NEKTAR_LIB_UTILITIES_FOUNDATIONS_POINTS_H
37 
38 #include <boost/core/ignore_unused.hpp>
39 
44 
45 //#include <LibUtilities/BasicUtils/BasicUtilsFwd.hpp> // for NekManager
48 
49 namespace Nektar
50 {
51 
52 namespace LibUtilities
53 {
54 // Need to add method to compute total number of points given dimension
55 // and number of points.
56 
57 /// Defines a specification for a set of points.
58 class PointsKey
59 {
60 public:
61  // Used for looking up the creator. The creator for number of points
62  // can generate for any number, so we want the same creator called
63  // for all number.
64  struct opLess
65  {
67  const PointsKey &rhs) const;
68  };
69 
70  /// Default constructor.
71  PointsKey(void)
73  m_factor(NekConstants::kNekUnsetDouble)
74  {
75  }
76 
77  /// Constructor defining the number and distribution of points.
78  PointsKey(const int &numpoints, const PointsType &pointstype,
80  : m_numpoints(numpoints), m_pointstype(pointstype), m_factor(factor)
81  {
82  }
83 
84  /// Destructor.
85  virtual ~PointsKey()
86  {
87  }
88 
89  /// Copy constructor.
90  PointsKey(const PointsKey &key)
91  {
92  *this = key; // defer to assignment operator
93  }
94 
96  {
99  m_factor = key.m_factor;
100 
101  return *this;
102  }
103 
104  inline unsigned int GetNumPoints() const
105  {
106  return m_numpoints;
107  }
108 
109  inline PointsType GetPointsType() const
110  {
111  return m_pointstype;
112  }
113 
114  inline NekDouble GetFactor() const
115  {
116  return m_factor;
117  }
118 
119  inline bool operator==(const PointsKey &key)
120  {
121 
122  if (fabs(m_factor - key.m_factor) < NekConstants::kNekZeroTol)
123  {
124  return (m_numpoints == key.m_numpoints &&
125  m_pointstype == key.m_pointstype);
126  }
127 
128  return false;
129  }
130 
131  inline bool operator==(const PointsKey *y)
132  {
133  return (*this == *y);
134  }
135 
136  inline bool operator!=(const PointsKey &y)
137  {
138  return (!(*this == y));
139  }
140 
141  inline bool operator!=(const PointsKey *y)
142  {
143  return (!(*this == *y));
144  }
145 
146  // If new points are added, this function must be modified
147  inline unsigned int GetPointsDim() const
148  {
149  int dimpoints = 1;
150 
151  switch (m_pointstype)
152  {
153  case eNodalTriElec:
154  case eNodalTriFekete:
156  case eNodalTriSPI:
157  case eNodalQuadElec:
158  dimpoints = 2;
159  break;
160 
161  case eNodalTetElec:
164  case eNodalPrismElec:
165  case eNodalHexElec:
166  dimpoints = 3;
167  break;
168 
169  default:
170  break;
171  }
172 
173  return dimpoints;
174  }
175 
176  // If new points are added, this function must be modified
177  inline unsigned int GetTotNumPoints() const
178  {
179  int totpoints = m_numpoints;
180 
181  switch (m_pointstype)
182  {
183  case eNodalTriElec:
184  case eNodalTriFekete:
186  totpoints = m_numpoints * (m_numpoints + 1) / 2;
187  break;
188  case eNodalTriSPI:
190  "This method cannot be implemented");
191  break;
192 
193  case eNodalQuadElec:
194  totpoints = m_numpoints * m_numpoints;
195  break;
196 
197  case eNodalTetElec:
199  totpoints =
200  m_numpoints * (m_numpoints + 1) * (m_numpoints + 2) / 6;
201  break;
202  case eNodalTetSPI:
204  "This method cannot be implemented");
205  break;
206 
208  case eNodalPrismElec:
209  totpoints = m_numpoints * m_numpoints * (m_numpoints + 1) / 2;
210  break;
211  case eNodalPrismSPI:
213  "This method cannot be implemented");
214  break;
215 
216  case eNodalHexElec:
217  totpoints = m_numpoints * m_numpoints * m_numpoints;
218  break;
219 
220  default:
221  break;
222  }
223 
224  return totpoints;
225  }
226 
227  LIB_UTILITIES_EXPORT friend bool operator==(const PointsKey &lhs,
228  const PointsKey &rhs);
229  LIB_UTILITIES_EXPORT friend bool operator<(const PointsKey &lhs,
230  const PointsKey &rhs);
232  const PointsKey &lhs, const PointsKey &rhs) const;
233 
234 protected:
235  unsigned int m_numpoints; //!< number of the points (as appropriately
236  //!< defined for PointsType)
237  PointsType m_pointstype; //!< Type of Points
238  NekDouble m_factor; //!< optional factor
239 private:
240 };
241 
243 
245  const PointsKey &rhs);
246 LIB_UTILITIES_EXPORT bool operator<(const PointsKey &lhs, const PointsKey &rhs);
247 LIB_UTILITIES_EXPORT std::ostream &operator<<(std::ostream &os,
248  const PointsKey &rhs);
249 
250 typedef std::vector<PointsKey> PointsKeyVector;
251 
252 /// Stores a set of points of datatype DataT, defined by a PointKey.
253 template <typename DataT> class Points
254 {
255 public:
256  typedef DataT DataType;
257  typedef std::shared_ptr<NekMatrix<DataType>> MatrixSharedPtrType;
258 
259  virtual ~Points()
260  {
261  }
262 
263  virtual void Initialize(void)
264  {
265  CalculatePoints();
269  }
270 
271  inline unsigned int GetPointsDim() const
272  {
273  return m_pointsKey.GetPointsDim();
274  }
275 
276  inline unsigned int GetNumPoints() const
277  {
278  return m_pointsKey.GetNumPoints();
279  }
280 
281  inline unsigned int GetTotNumPoints() const
282  {
283  return m_pointsKey.GetTotNumPoints();
284  }
285 
286  inline PointsType GetPointsType() const
287  {
288  return m_pointsKey.GetPointsType();
289  }
290 
291  inline const Array<OneD, const DataType> &GetZ() const
292  {
293  return m_points[0];
294  }
295 
296  inline const Array<OneD, const DataType> &GetW() const
297  {
298  return m_weights;
299  }
300 
303  {
304  z = m_points[0];
305  w = m_weights;
306  }
307 
309  {
310  return m_bcweights;
311  }
312 
314  {
315  x = m_points[0];
316  }
317 
320  {
321  x = m_points[0];
322  y = m_points[1];
323  }
324 
328  {
329  x = m_points[0];
330  y = m_points[1];
331  z = m_points[2];
332  }
333 
334  inline const MatrixSharedPtrType &GetD(Direction dir = xDir) const
335  {
336  return m_derivmatrix[(int)dir];
337  }
338 
339  virtual const MatrixSharedPtrType GetI(const PointsKey &key)
340  {
341  boost::ignore_unused(key);
342  NEKERROR(ErrorUtil::efatal, "Method not implemented ");
343  std::shared_ptr<NekMatrix<NekDouble>> returnval(
344  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr());
345  return returnval;
346  }
347 
349  {
350  boost::ignore_unused(x);
351  NEKERROR(ErrorUtil::efatal, "Method not implemented");
352  std::shared_ptr<NekMatrix<NekDouble>> returnval(
353  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr());
354  return returnval;
355  }
356 
357  virtual const MatrixSharedPtrType GetI(unsigned int,
359  {
360  boost::ignore_unused(x);
361  NEKERROR(ErrorUtil::efatal, "Method not implemented");
362  std::shared_ptr<NekMatrix<NekDouble>> returnval(
363  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr());
364  return returnval;
365  }
366 
369  {
370  boost::ignore_unused(x, y);
371  NEKERROR(ErrorUtil::efatal, "Method not implemented");
372  std::shared_ptr<NekMatrix<NekDouble>> returnval(
373  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr());
374  return returnval;
375  }
376 
380  {
381  boost::ignore_unused(x, y, z);
382  NEKERROR(ErrorUtil::efatal, "Method not implemented");
383  std::shared_ptr<NekMatrix<NekDouble>> returnval(
384  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr());
385  return returnval;
386  }
387 
389  const PointsKey &pkey)
390  {
391  boost::ignore_unused(pkey);
392  NEKERROR(ErrorUtil::efatal, "Method not implemented ");
393  std::shared_ptr<NekMatrix<NekDouble>> returnval(
394  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr());
395  return returnval;
396  }
397 
398 protected:
399  /// Points type for this points distributions.
401  /// Storage for the point locations, allowing for up to a 3D points
402  /// storage.
404  /// Quadrature weights for the weights.
406  /// Barycentric weights.
408  /// Derivative matrices.
414 
415  virtual void CalculatePoints()
416  {
417  unsigned int pointsDim = GetPointsDim();
418  unsigned int totNumPoints = GetTotNumPoints();
419 
420  for (unsigned int i = 0; i < pointsDim; ++i)
421  {
422  m_points[i] = Array<OneD, DataType>(totNumPoints);
423  }
424  }
425 
426  virtual void CalculateWeights()
427  {
429  }
430 
431  /**
432  * @brief This function calculates the barycentric weights used for
433  * enhanced interpolation speed.
434  *
435  * For the points distribution \f$ z_i \f$ with \f% 1\leq z_i \leq N
436  * \f$, the barycentric weights are computed as:
437  *
438  * \f[
439  * b_i=\prod_{\substack{1\leq j\leq N\\ i\neq j}} \frac{1}{z_i-z_j}
440  * \f]
441  */
442  virtual void CalculateBaryWeights()
443  {
444  const unsigned int totNumPoints = m_pointsKey.GetNumPoints();
445  m_bcweights = Array<OneD, DataType>(totNumPoints, 1.0);
446 
448 
449  for (unsigned int i = 0; i < totNumPoints; ++i)
450  {
451  for (unsigned int j = 0; j < totNumPoints; ++j)
452  {
453  if (i == j)
454  {
455  continue;
456  }
457 
458  m_bcweights[i] *= (z[i] - z[j]);
459  }
460 
461  m_bcweights[i] = 1.0 / m_bcweights[i];
462  }
463  }
464 
465  virtual void CalculateDerivMatrix()
466  {
467  int totNumPoints = GetTotNumPoints();
468  for (unsigned int i = 0; i < m_pointsKey.GetPointsDim(); ++i)
469  {
470  m_derivmatrix[i] =
471  MemoryManager<NekMatrix<DataType>>::AllocateSharedPtr(
472  totNumPoints, totNumPoints);
473  }
474  }
475 
476  Points(const PointsKey &key) : m_pointsKey(key)
477  {
478  }
479 
480 private:
481  // These should never be called
482  Points(const Points &pts)
483  {
484  boost::ignore_unused(pts);
486  "Copy Constructor for Points should not be called");
487  }
489  {
491  "Default Constructor for Points should not be called");
492  }
493 };
494 
495 } // namespace LibUtilities
496 } // namespace Nektar
497 
498 #endif // NEKTAR_LIB_UTILITIES_FOUNDATIONS_POINTS_H
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define LIB_UTILITIES_EXPORT
1D Array of constant elements with garbage collection and bounds checking.
Definition: SharedArray.hpp:58
Stores a set of points of datatype DataT, defined by a PointKey.
Definition: Points.h:254
Array< OneD, DataType > m_points[3]
Storage for the point locations, allowing for up to a 3D points storage.
Definition: Points.h:403
MatrixSharedPtrType m_derivmatrix[3]
Derivative matrices.
Definition: Points.h:409
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess > m_InterpManager
Definition: Points.h:411
const MatrixSharedPtrType & GetD(Direction dir=xDir) const
Definition: Points.h:334
virtual void CalculateWeights()
Definition: Points.h:426
void GetPoints(Array< OneD, const DataType > &x, Array< OneD, const DataType > &y) const
Definition: Points.h:318
virtual const MatrixSharedPtrType GetGalerkinProjection(const PointsKey &pkey)
Definition: Points.h:388
virtual void CalculateDerivMatrix()
Definition: Points.h:465
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess > m_GalerkinProjectionManager
Definition: Points.h:413
PointsKey m_pointsKey
Points type for this points distributions.
Definition: Points.h:400
void GetZW(Array< OneD, const DataType > &z, Array< OneD, const DataType > &w) const
Definition: Points.h:301
unsigned int GetNumPoints() const
Definition: Points.h:276
void GetPoints(Array< OneD, const DataType > &x) const
Definition: Points.h:313
void GetPoints(Array< OneD, const DataType > &x, Array< OneD, const DataType > &y, Array< OneD, const DataType > &z) const
Definition: Points.h:325
PointsType GetPointsType() const
Definition: Points.h:286
Points(const PointsKey &key)
Definition: Points.h:476
std::shared_ptr< NekMatrix< DataType > > MatrixSharedPtrType
Definition: Points.h:257
virtual const MatrixSharedPtrType GetI(const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y, const Array< OneD, const DataType > &z)
Definition: Points.h:377
virtual void Initialize(void)
Definition: Points.h:263
const Array< OneD, const DataType > & GetW() const
Definition: Points.h:296
virtual void CalculatePoints()
Definition: Points.h:415
virtual const MatrixSharedPtrType GetI(const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y)
Definition: Points.h:367
virtual const MatrixSharedPtrType GetI(const Array< OneD, const DataType > &x)
Definition: Points.h:348
unsigned int GetTotNumPoints() const
Definition: Points.h:281
unsigned int GetPointsDim() const
Definition: Points.h:271
Array< OneD, DataType > m_weights
Quadrature weights for the weights.
Definition: Points.h:405
const Array< OneD, const DataType > & GetZ() const
Definition: Points.h:291
virtual void CalculateBaryWeights()
This function calculates the barycentric weights used for enhanced interpolation speed.
Definition: Points.h:442
Points(const Points &pts)
Definition: Points.h:482
Array< OneD, DataType > m_bcweights
Barycentric weights.
Definition: Points.h:407
virtual const MatrixSharedPtrType GetI(unsigned int, const Array< OneD, const DataType > &x)
Definition: Points.h:357
virtual const MatrixSharedPtrType GetI(const PointsKey &key)
Definition: Points.h:339
const Array< OneD, const NekDouble > & GetBaryWeights() const
Definition: Points.h:308
Defines a specification for a set of points.
Definition: Points.h:59
unsigned int m_numpoints
number of the points (as appropriately defined for PointsType)
Definition: Points.h:235
unsigned int GetTotNumPoints() const
Definition: Points.h:177
bool operator==(const PointsKey *y)
Definition: Points.h:131
PointsType GetPointsType() const
Definition: Points.h:109
friend bool operator<(const PointsKey &lhs, const PointsKey &rhs)
PointsType m_pointstype
Type of Points.
Definition: Points.h:237
PointsKey & operator=(const PointsKey &key)
Definition: Points.h:95
PointsKey(const int &numpoints, const PointsType &pointstype, const NekDouble factor=NekConstants::kNekUnsetDouble)
Constructor defining the number and distribution of points.
Definition: Points.h:78
unsigned int GetPointsDim() const
Definition: Points.h:147
bool operator==(const PointsKey &key)
Definition: Points.h:119
PointsKey(const PointsKey &key)
Copy constructor.
Definition: Points.h:90
NekDouble GetFactor() const
Definition: Points.h:114
NekDouble m_factor
optional factor
Definition: Points.h:238
unsigned int GetNumPoints() const
Definition: Points.h:104
virtual ~PointsKey()
Destructor.
Definition: Points.h:85
bool operator!=(const PointsKey &y)
Definition: Points.h:136
PointsKey(void)
Default constructor.
Definition: Points.h:71
bool operator!=(const PointsKey *y)
Definition: Points.h:141
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
bool operator==(const BasisKey &x, const BasisKey &y)
bool operator<(const BasisKey &lhs, const BasisKey &rhs)
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
@ eNodalPrismEvenlySpaced
3D Evenly-spaced points on a Prism
Definition: PointsType.h:88
@ eNodalTriFekete
2D Nodal Fekete Points on a Triangle
Definition: PointsType.h:84
@ eNodalPrismSPI
3D prism SPI
Definition: PointsType.h:94
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:83
@ eNodalTriEvenlySpaced
2D Evenly-spaced points on a Triangle
Definition: PointsType.h:85
@ eNodalHexElec
3D GLL for hex
Definition: PointsType.h:96
@ eNodalQuadElec
2D GLL for quad
Definition: PointsType.h:95
@ eNodalTetEvenlySpaced
3D Evenly-spaced points on a Tetrahedron
Definition: PointsType.h:86
@ eNodalTetSPI
3D Nodal Symmetric positive internal tet (Whitherden, Vincent)
Definition: PointsType.h:92
@ eNodalPrismElec
3D electrostatically spaced points on a Prism
Definition: PointsType.h:89
@ eNodalTetElec
3D Nodal Electrostatic Points on a Tetrahedron
Definition: PointsType.h:87
@ eNodalTriSPI
2D Nodal Symmetric positive internal triangle (Whitherden, Vincent)
Definition: PointsType.h:90
static const PointsKey NullPointsKey(0, eNoPointsType)
static const NekDouble kNekUnsetDouble
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
bool operator()(const PointsKey &lhs, const PointsKey &rhs) const