Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
NodalUtil.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalUtil.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: 2D Nodal Triangle Fekete Utilities --
33 // Basis function, Interpolation, Integral, Derivation, etc
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iomanip>
38 #include <limits>
39 
42 
47 
48 namespace Nektar
49 {
50 namespace LibUtilities
51 {
52 
53 /**
54  * @brief Obtain the integration weights for the given nodal distribution.
55  *
56  * This routine constructs the integration weights for the given nodal
57  * distribution inside NodalUtil::m_xi. To do this we solve the linear system
58  * \f$ \mathbf{V}^\top \mathbf{w} = \mathbf{g} \f$ where \f$ \mathbf{g}_i =
59  * \int_\Omega \psi_i(x)\, dx \f$, and we use the fact that under the definition
60  * of the orthogonal basis for each element type, \f$ \mathbf{g}_i = 0 \f$ for
61  * \f$ i > 0 \f$. We use NodalUtil::v_ModeZeroIntegral to return the analytic
62  * value of \f$ \mathbf{g}_0 \f$.
63  *
64  * @return Vector of integration weights for the integration points.
65  */
67 {
68  // Get number of modes in orthogonal basis
69  int numModes = v_NumModes();
70 
71  // If we have the same number of nodes as modes, then we can solve the
72  // linear system V^T w = (1,0,...)
73  if (numModes == m_numPoints)
74  {
76  g(0) = v_ModeZeroIntegral();
77 
79 
80  // Solve the system V^T w = g to obtain weights.
81  LinearSystem matL(V);
82  return matL.SolveTranspose(g);
83  }
84  else
85  {
86  // System is either over- or under-determined. Need to do least squares
87  // here using SVD.
88  return NekVector<NekDouble>();
89  }
90 }
91 
92 /**
93  * @brief Return the Vandermonde matrix for the nodal distribution.
94  *
95  * This routine constructs and returns the Vandermonde matrix \f$\mathbf{V}\f$,
96  * with each entry as \f$\mathbf{V}_{ij} = (\psi_i(\xi_j))\f$ where \f$ \psi_i
97  * \f$ is the orthogonal basis obtained through the abstract function
98  * NodalUtil::v_OrthoBasis.
99  *
100  * @return The Vandermonde matrix.
101  */
103 {
104  int rows = m_numPoints, cols = v_NumModes();
105  SharedMatrix matV = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
106  rows, cols, 0.0);
107 
108  for (int j = 0; j < cols; ++j)
109  {
111  for (int i = 0; i < rows; ++i)
112  {
113  (*matV)(i, j) = col[i];
114  }
115  }
116 
117  return matV;
118 }
119 /**
120  * @brief Return the Vandermonde matrix of the derivative of the basis functions
121  * for the nodal distribution.
122  *
123  * This routine constructs and returns the Vandermonde matrix for the derivative
124  * of the basis functions \f$\mathbf{V}_d\f$ for coordinate directions \f$ 0
125  * \leq d \leq 2 \f$, with each entry as \f$\mathbf{V}_{ij} = (\partial_d
126  * \psi_i(\xi_j))\f$ where \f$ \partial_d\psi_i \f$ is the derivative of the
127  * orthogonal basis obtained through the abstract function
128  * NodalUtil::v_OrthoBasisDeriv.
129  *
130  * @param dir Direction of derivative in the standard element.
131  *
132  * @return Vandermonde matrix corresponding with derivative of the basis
133  * functions in direction @p dir.
134  */
136 {
137  int rows = m_numPoints, cols = v_NumModes();
138  SharedMatrix matV = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
139  rows, cols, 0.0);
140 
141  for (int j = 0; j < cols; ++j)
142  {
144  for (int i = 0; i < rows; ++i)
145  {
146  (*matV)(i, j) = col[i];
147  }
148  }
149 
150  return matV;
151 }
152 
153 /**
154  * @brief Return the derivative matrix for the nodal distribution.
155  *
156  * This routine constructs and returns the derivative matrices
157  * \f$\mathbf{D}_d\f$ for coordinate directions \f$ 0 \leq d \leq 2 \f$, which
158  * can be used to evaluate the derivative of a nodal expansion at the points
159  * defined by NodalUtil::m_xi. These are calculated as \f$ \mathbf{D}_d =
160  * \mathbf{V}_d \mathbf{V}^{-1} \f$, where \f$ \mathbf{V}_d \f$ is the
161  * derivative Vandermonde matrix and \f$ \mathbf{V} \f$ is the Vandermonde
162  * matrix.
163  *
164  * @param dir Coordinate direction in which to evaluate the derivative.
165  *
166  * @return The derivative matrix for direction @p dir.
167  */
169 {
172  SharedMatrix D = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
173  V->GetRows(), V->GetColumns(), 0.0);
174 
175  V->Invert();
176 
177  *D = (*Vd) * (*V);
178 
179  return D;
180 }
181 
182 /**
183  * @brief Construct the interpolation matrix used to evaluate the basis at the
184  * points @p xi inside the element.
185  *
186  * This routine returns a matrix \f$ \mathbf{I}(\mathbf{a}) \f$ that can be used
187  * to evaluate the nodal basis at the points defined by the parameter @p xi,
188  * which is denoted by \f$ \mathbf{a} = (a_1, \dots, a_N) \f$ and \f$ N \f$ is
189  * the number of points in @p xi.
190  *
191  * In particular, if the array \f$ \mathbf{u} \f$ with components \f$
192  * \mathbf{u}_i = u^\delta(\xi_i) \f$ represents the polynomial approximation of
193  * a function \f$ u \f$ evaluated at the nodal points NodalUtil::m_xi, then the
194  * evaluation of \f$ u^\delta \f$ evaluated at the input points \f$ \mathbf{a}
195  * \f$ is given by \f$ \mathbf{I}(\mathbf{a})\mathbf{u} \f$.
196  *
197  * @param xi An array of first size number of spatial dimensions \f$ d \f$ and
198  * secondary size the number of points to interpolate.
199  *
200  * @return The interpolation matrix for the points @p xi.
201  */
204 {
205  boost::shared_ptr<NodalUtil> subUtil = v_CreateUtil(xi);
206  SharedMatrix matS = GetVandermonde();
207  SharedMatrix matT = subUtil->GetVandermonde();
208  SharedMatrix D = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
209  matT->GetRows(), matS->GetColumns(), 0.0);
210 
211  matS->Invert();
212 
213  *D = (*matT) * (*matS);
214  return D;
215 }
216 
217 /**
218  * @brief Construct the nodal utility class for a triangle.
219  *
220  * The constructor of this class sets up two member variables used in the
221  * evaluation of the orthogonal basis:
222  *
223  * - NodalUtilTriangle::m_eta is used to construct the collapsed coordinate
224  * locations of the nodal points \f$ (\eta_1, \eta_2) \f$ inside the square
225  * \f$[-1,1]^2\f$ on which the orthogonal basis functions are defined.
226  * - NodalUtilTriangle::m_ordering constructs a mapping from the index set \f$ I
227  * = \{ (i,j)\ |\ 0\leq i,j \leq P, i+j \leq P \}\f$ to an ordering \f$ 0 \leq
228  * m(ij) \leq (P+1)(P+2)/2 \f$ that defines the monomials \f$ \xi_1^i \xi_2^j
229  * \f$ that span the triangular space. This is then used to calculate which
230  * \f$ (i,j) \f$ pair corresponding to a column of the Vandermonde matrix when
231  * calculating the orthogonal polynomials.
232  *
233  * @param degree Polynomial order of this nodal triangle.
234  * @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
235  * element.
236  * @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
237  * element.
238  */
242  : NodalUtil(degree, 2), m_eta(2)
243 {
244  // Set up parent variables.
245  m_numPoints = r.num_elements();
246  m_xi[0] = r;
247  m_xi[1] = s;
248 
249  // Construct a mapping (i,j) -> m from the triangular tensor product space
250  // (i,j) to a single ordering m.
251  for (int i = 0; i <= m_degree; ++i)
252  {
253  for (int j = 0; j <= m_degree - i; ++j)
254  {
255  m_ordering.push_back(std::make_pair(i,j));
256  }
257  }
258 
259  // Calculate collapsed coordinates from r/s values
262 
263  for (int i = 0; i < m_numPoints; ++i)
264  {
265  if (fabs(m_xi[1][i]-1.0) < NekConstants::kNekZeroTol)
266  {
267  m_eta[0][i] = -1.0;
268  m_eta[1][i] = 1.0;
269  }
270  else
271  {
272  m_eta[0][i] = 2*(1+m_xi[0][i])/(1-m_xi[1][i])-1.0;
273  m_eta[1][i] = m_xi[1][i];
274  }
275  }
276 }
277 
278 /**
279  * @brief Return the value of the modal functions for the triangular element at
280  * the nodal points #m_xi for a given mode.
281  *
282  * In a triangle, we use the orthogonal basis
283  *
284  * \f[
285  * \psi_{m(ij)} = \sqrt{2} P^{(0,0)}_i(\xi_1) P_j^{(2i+1,0)}(\xi_2) (1-\xi_2)^i
286  * \f]
287  *
288  * where \f$ m(ij) \f$ is the mapping defined in NodalUtilTriangle::m_ordering
289  * and \f$ J_n^{(\alpha,\beta)}(z) \f$ denotes the standard Jacobi polynomial.
290  *
291  * @param mode The mode of the orthogonal basis to evaluate.
292  *
293  * @return Vector containing orthogonal basis evaluated at the points #m_xi.
294  */
296 {
297  std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
298  std::pair<int, int> modes = m_ordering[mode];
299 
300  // Calculate Jacobi polynomials
302  m_numPoints, &m_eta[0][0], &jacobi_i[0], NULL, modes.first, 0.0, 0.0);
304  m_numPoints, &m_eta[1][0], &jacobi_j[0], NULL, modes.second,
305  2.0 * modes.first + 1.0, 0.0);
306 
308  NekDouble sqrt2 = sqrt(2.0);
309 
310  for (int i = 0; i < m_numPoints; ++i)
311  {
312  ret[i] = sqrt2 * jacobi_i[i] * jacobi_j[i] *
313  pow(1.0 - m_eta[1][i], modes.first);
314  }
315 
316  return ret;
317 }
318 
319 /**
320  * @brief Return the value of the derivative of the modal functions for the
321  * triangular element at the nodal points #m_xi for a given mode.
322  *
323  * Note that this routine must use the chain rule combined with the collapsed
324  * coordinate derivatives as described in Sherwin & Karniadakis (2nd edition),
325  * pg 150.
326  *
327  * @param dir Coordinate direction in which to evaluate the derivative.
328  * @param mode The mode of the orthogonal basis to evaluate.
329  *
330  * @return Vector containing the derivative of the orthogonal basis evaluated at
331  * the points #m_xi.
332  */
334  const int dir, const int mode)
335 {
336  std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
337  std::vector<NekDouble> jacobi_di(m_numPoints), jacobi_dj(m_numPoints);
338  std::pair<int, int> modes = m_ordering[mode];
339 
340  // Calculate Jacobi polynomials and their derivatives. Note that we use both
341  // jacobfd and jacobd since jacobfd is only valid for derivatives in the
342  // open interval (-1,1).
344  m_numPoints, &m_eta[0][0], &jacobi_i[0], NULL, modes.first, 0.0,
345  0.0);
347  m_numPoints, &m_eta[1][0], &jacobi_j[0], NULL, modes.second,
348  2.0*modes.first + 1.0, 0.0);
350  m_numPoints, &m_eta[0][0], &jacobi_di[0], modes.first, 0.0, 0.0);
352  m_numPoints, &m_eta[1][0], &jacobi_dj[0], modes.second,
353  2.0*modes.first + 1.0, 0.0);
354 
356  NekDouble sqrt2 = sqrt(2.0);
357 
358  if (dir == 0)
359  {
360  // d/d(\xi_1) = 2/(1-\eta_2) d/d(\eta_1)
361  for (int i = 0; i < m_numPoints; ++i)
362  {
363  ret[i] = 2.0 * sqrt2 * jacobi_di[i] * jacobi_j[i];
364  if (modes.first > 0)
365  {
366  ret[i] *= pow(1.0 - m_eta[1][i], modes.first - 1.0);
367  }
368  }
369  }
370  else
371  {
372  // d/d(\xi_2) = 2(1+\eta_1)/(1-\eta_2) d/d(\eta_1) + d/d(eta_2)
373  for (int i = 0; i < m_numPoints; ++i)
374  {
375  ret[i] = (1 + m_eta[0][i]) * sqrt2 * jacobi_di[i] * jacobi_j[i];
376  if (modes.first > 0)
377  {
378  ret[i] *= pow(1.0 - m_eta[1][i], modes.first - 1.0);
379  }
380 
381  NekDouble tmp = jacobi_dj[i] * pow(1.0 - m_eta[1][i], modes.first);
382  if (modes.first > 0)
383  {
384  tmp -= modes.first * jacobi_j[i] *
385  pow(1.0 - m_eta[1][i], modes.first-1);
386  }
387 
388  ret[i] += sqrt2 * jacobi_i[i] * tmp;
389  }
390  }
391 
392  return ret;
393 }
394 
395 /**
396  * @brief Construct the nodal utility class for a tetrahedron.
397  *
398  * The constructor of this class sets up two member variables used in the
399  * evaluation of the orthogonal basis:
400  *
401  * - NodalUtilTetrahedron::m_eta is used to construct the collapsed coordinate
402  * locations of the nodal points \f$ (\eta_1, \eta_2, \eta_3) \f$ inside the
403  * cube \f$[-1,1]^3\f$ on which the orthogonal basis functions are defined.
404  * - NodalUtilTetrahedron::m_ordering constructs a mapping from the index set
405  * \f$ I = \{ (i,j,k)\ |\ 0\leq i,j,k \leq P, i+j \leq P, i+j+k \leq P \}\f$
406  * to an ordering \f$ 0 \leq m(ijk) \leq (P+1)(P+2)(P+3)/6 \f$ that defines
407  * the monomials \f$ \xi_1^i \xi_2^j \xi_3^k \f$ that span the tetrahedral
408  * space. This is then used to calculate which \f$ (i,j,k) \f$ triple
409  * (represented as a boost tuple) corresponding to a column of the Vandermonde
410  * matrix when calculating the orthogonal polynomials.
411  *
412  * @param degree Polynomial order of this nodal tetrahedron
413  * @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
414  * element.
415  * @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
416  * element.
417  * @param t \f$ \xi_3 \f$-coordinates of nodal points in the standard
418  * element.
419  */
424  : NodalUtil(degree, 3), m_eta(3)
425 {
426  m_numPoints = r.num_elements();
427  m_xi[0] = r;
428  m_xi[1] = s;
429  m_xi[2] = t;
430 
431  for (int i = 0; i <= m_degree; ++i)
432  {
433  for (int j = 0; j <= m_degree - i; ++j)
434  {
435  for (int k = 0; k <= m_degree - i - j; ++k)
436  {
437  m_ordering.push_back(Mode(i, j, k));
438  }
439  }
440  }
441 
442  // Calculate collapsed coordinates from r/s values
446 
447  for (int i = 0; i < m_numPoints; ++i)
448  {
449  if (fabs(m_xi[2][i] - 1.0) < NekConstants::kNekZeroTol)
450  {
451  // Very top point of the tetrahedron
452  m_eta[0][i] = -1.0;
453  m_eta[1][i] = -1.0;
454  m_eta[2][i] = m_xi[2][i];
455  }
456  else
457  {
458  if (fabs(m_xi[1][i] - 1.0) < NekConstants::kNekZeroTol)
459  {
460  // Distant diagonal edge shared by all eta_x coordinate planes:
461  // the xi_y == -xi_z line
462  m_eta[0][i] = -1.0;
463  }
464  else if (fabs(m_xi[1][i] + m_xi[2][i]) < NekConstants::kNekZeroTol)
465  {
466  m_eta[0][i] = -1.0;
467  }
468  else
469  {
470  m_eta[0][i] = 2.0 * (1.0 + m_xi[0][i]) /
471  (-m_xi[1][i] - m_xi[2][i]) - 1.0;
472  }
473  m_eta[1][i] = 2.0 * (1.0 + m_xi[1][i]) / (1.0 - m_xi[2][i]) - 1.0;
474  m_eta[2][i] = m_xi[2][i];
475  }
476  }
477 }
478 
479 /**
480  * @brief Return the value of the modal functions for the tetrahedral element at
481  * the nodal points #m_xi for a given mode.
482  *
483  * In a tetrahedron, we use the orthogonal basis
484  *
485  * \f[ \psi_{m(ijk)} = \sqrt{8} P^{(0,0)}_i(\xi_1) P_j^{(2i+1,0)}(\xi_2)
486  * P_k^{(2i+2j+2,0)}(\xi_3) (1-\xi_2)^i (1-\xi_3)^{i+j} \f]
487  *
488  * where \f$ m(ijk) \f$ is the mapping defined in #m_ordering and \f$
489  * J_n^{(\alpha,\beta)}(z) \f$ denotes the standard Jacobi polynomial.
490  *
491  * @param mode The mode of the orthogonal basis to evaluate.
492  *
493  * @return Vector containing orthogonal basis evaluated at the points #m_xi.
494  */
496 {
497  std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
498  std::vector<NekDouble> jacC(m_numPoints);
499  Mode modes = m_ordering[mode];
500 
501  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
502 
503  // Calculate Jacobi polynomials
505  m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
507  m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 2.0 * I + 1.0, 0.0);
509  m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * (I+J) + 2.0, 0.0);
510 
512  NekDouble sqrt8 = sqrt(8.0);
513 
514  for (int i = 0; i < m_numPoints; ++i)
515  {
516  ret[i] = sqrt8 * jacA[i] * jacB[i] * jacC[i] *
517  pow(1.0 - m_eta[1][i], I) * pow(1.0 - m_eta[2][i], I + J);
518  }
519 
520  return ret;
521 }
522 
523 /**
524  * @brief Return the value of the derivative of the modal functions for the
525  * tetrahedral element at the nodal points #m_xi for a given mode.
526  *
527  * Note that this routine must use the chain rule combined with the collapsed
528  * coordinate derivatives as described in Sherwin & Karniadakis (2nd edition),
529  * pg 152.
530  *
531  * @param dir Coordinate direction in which to evaluate the derivative.
532  * @param mode The mode of the orthogonal basis to evaluate.
533  *
534  * @return Vector containing the derivative of the orthogonal basis evaluated at
535  * the points #m_xi.
536  */
538  const int dir, const int mode)
539 {
540  std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
541  std::vector<NekDouble> jacC(m_numPoints);
542  std::vector<NekDouble> jacDerivA(m_numPoints), jacDerivB(m_numPoints);
543  std::vector<NekDouble> jacDerivC(m_numPoints);
544  Mode modes = m_ordering[mode];
545 
546  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
547 
548  // Calculate Jacobi polynomials
550  m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
552  m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 2.0 * I + 1.0, 0.0);
554  m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * (I+J) + 2.0, 0.0);
556  m_numPoints, &m_eta[0][0], &jacDerivA[0], I, 0.0, 0.0);
558  m_numPoints, &m_eta[1][0], &jacDerivB[0], J, 2.0 * I + 1.0, 0.0);
560  m_numPoints, &m_eta[2][0], &jacDerivC[0], K, 2.0 * (I+J) + 2.0, 0.0);
561 
563  NekDouble sqrt8 = sqrt(8.0);
564 
565  // Always compute x-derivative since this term appears in the latter two
566  // terms.
567  for (int i = 0; i < m_numPoints; ++i)
568  {
569  ret[i] = 4.0 * sqrt8 * jacDerivA[i] * jacB[i] * jacC[i];
570 
571  if (I > 0)
572  {
573  ret[i] *= pow(1 - m_eta[1][i], I - 1);
574  }
575 
576  if (I + J > 0)
577  {
578  ret[i] *= pow(1 - m_eta[2][i], I + J - 1);
579  }
580  }
581 
582  if (dir >= 1)
583  {
584  // Multiply by (1+a)/2
585  NekVector<NekDouble> tmp(m_numPoints);
586 
587  for (int i = 0; i < m_numPoints; ++i)
588  {
589  ret[i] *= 0.5 * (m_eta[0][i] + 1.0);
590 
591  tmp[i] = 2.0 * sqrt8 * jacA[i] * jacC[i];
592  if (I + J > 0)
593  {
594  tmp[i] *= pow(1.0 - m_eta[2][i], I + J - 1);
595  }
596 
597  NekDouble tmp2 = jacDerivB[i] * pow(1.0 - m_eta[1][i], I);
598  if (I > 0)
599  {
600  tmp2 -= I * jacB[i] * pow(1.0 - m_eta[1][i], I - 1);
601  }
602 
603  tmp[i] *= tmp2;
604  }
605 
606  if (dir == 1)
607  {
608  ret += tmp;
609  return ret;
610  }
611 
612  for (int i = 0; i < m_numPoints; ++i)
613  {
614  ret[i] += 0.5 * (1.0 + m_eta[1][i]) * tmp[i];
615 
616  NekDouble tmp2 = jacDerivC[i] * pow(1.0 - m_eta[2][i], I + J);
617  if (I + J > 0)
618  {
619  tmp2 -= jacC[i] * (I + J) * pow(1.0 - m_eta[2][i], I + J - 1);
620  }
621 
622  ret[i] += sqrt8 * jacA[i] * jacB[i] * pow(1.0 - m_eta[1][i], I) *
623  tmp2;
624  }
625  }
626 
627  return ret;
628 }
629 
630 /**
631  * @brief Construct the nodal utility class for a prism.
632  *
633  * The constructor of this class sets up two member variables used in the
634  * evaluation of the orthogonal basis:
635  *
636  * - NodalUtilPrism::m_eta is used to construct the collapsed coordinate
637  * locations of the nodal points \f$ (\eta_1, \eta_2, \eta_3) \f$ inside the
638  * cube \f$[-1,1]^3\f$ on which the orthogonal basis functions are defined.
639  * - NodalUtilPrism::m_ordering constructs a mapping from the index set
640  * \f$ I = \{ (i,j,k)\ |\ 0\leq i,j,k \leq P, i+k \leq P \}\f$ to an ordering
641  * \f$ 0 \leq m(ijk) \leq (P+1)(P+1)(P+2)/2 \f$ that defines the monomials \f$
642  * \xi_1^i \xi_2^j \xi_3^k \f$ that span the prismatic space. This is then
643  * used to calculate which \f$ (i,j,k) \f$ triple (represented as a boost
644  * tuple) corresponding to a column of the Vandermonde matrix when calculating
645  * the orthogonal polynomials.
646  *
647  * @param degree Polynomial order of this nodal tetrahedron
648  * @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
649  * element.
650  * @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
651  * element.
652  * @param t \f$ \xi_3 \f$-coordinates of nodal points in the standard
653  * element.
654  */
659  : NodalUtil(degree, 3), m_eta(3)
660 {
661  m_numPoints = r.num_elements();
662  m_xi[0] = r;
663  m_xi[1] = s;
664  m_xi[2] = t;
665 
666  for (int i = 0; i <= m_degree; ++i)
667  {
668  for (int j = 0; j <= m_degree; ++j)
669  {
670  for (int k = 0; k <= m_degree - i; ++k)
671  {
672  m_ordering.push_back(Mode(i, j, k));
673  }
674  }
675  }
676 
677  // Calculate collapsed coordinates from r/s values
681 
682  for (int i = 0; i < m_numPoints; ++i)
683  {
684  if (fabs(m_xi[2][i] - 1.0) < NekConstants::kNekZeroTol)
685  {
686  // Very top point of the prism
687  m_eta[0][i] = -1.0;
688  m_eta[1][i] = m_xi[1][i];
689  m_eta[2][i] = 1.0;
690  }
691  else
692  {
693  // Third basis function collapsed to "pr" direction instead of "qr"
694  // direction
695  m_eta[0][i] = 2.0*(1.0 + m_xi[0][i])/(1.0 - m_xi[2][i]) - 1.0;
696  m_eta[1][i] = m_xi[1][i];
697  m_eta[2][i] = m_xi[2][i];
698  }
699  }
700 }
701 
702 /**
703  * @brief Return the value of the modal functions for the prismatic element at
704  * the nodal points #m_xi for a given mode.
705  *
706  * In a prism, we use the orthogonal basis
707  *
708  * \f[ \psi_{m(ijk)} = \sqrt{2} P^{(0,0)}_i(\xi_1) P_j^{(0,0)}(\xi_2)
709  * P_k^{(2i+1,0)}(\xi_3) (1-\xi_3)^i \f]
710  *
711  * where \f$ m(ijk) \f$ is the mapping defined in #m_ordering and \f$
712  * J_n^{(\alpha,\beta)}(z) \f$ denotes the standard Jacobi polynomial.
713  *
714  * @param mode The mode of the orthogonal basis to evaluate.
715  *
716  * @return Vector containing orthogonal basis evaluated at the points #m_xi.
717  */
719 {
720  std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
721  std::vector<NekDouble> jacC(m_numPoints);
722  Mode modes = m_ordering[mode];
723 
724  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
725 
726  // Calculate Jacobi polynomials
728  m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
730  m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 0.0, 0.0);
732  m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * I + 1.0, 0.0);
733 
735  NekDouble sqrt2 = sqrt(2.0);
736 
737  for (int i = 0; i < m_numPoints; ++i)
738  {
739  ret[i] = sqrt2 * jacA[i] * jacB[i] * jacC[i] *
740  pow(1.0 - m_eta[2][i], I);
741  }
742 
743  return ret;
744 }
745 
746 /**
747  * @brief Return the value of the derivative of the modal functions for the
748  * prismatic element at the nodal points #m_xi for a given mode.
749  *
750  * Note that this routine must use the chain rule combined with the collapsed
751  * coordinate derivatives as described in Sherwin & Karniadakis (2nd edition),
752  * pg 152.
753  *
754  * @param mode The mode of the orthogonal basis to evaluate.
755  * @param dir Coordinate direction in which to evaluate the derivative.
756  *
757  * @return Vector containing the derivative of the orthogonal basis evaluated at
758  * the points #m_xi.
759  */
761  const int dir, const int mode)
762 {
763  std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
764  std::vector<NekDouble> jacC(m_numPoints);
765  std::vector<NekDouble> jacDerivA(m_numPoints), jacDerivB(m_numPoints);
766  std::vector<NekDouble> jacDerivC(m_numPoints);
767  Mode modes = m_ordering[mode];
768 
769  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
770 
771  // Calculate Jacobi polynomials
773  m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
775  m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 0.0, 0.0);
777  m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * I + 1.0, 0.0);
779  m_numPoints, &m_eta[0][0], &jacDerivA[0], I, 0.0, 0.0);
781  m_numPoints, &m_eta[1][0], &jacDerivB[0], J, 0.0, 0.0);
783  m_numPoints, &m_eta[2][0], &jacDerivC[0], K, 2.0 * I + 1.0, 0.0);
784 
786  NekDouble sqrt2 = sqrt(2.0);
787 
788  if (dir == 1)
789  {
790  for (int i = 0; i < m_numPoints; ++i)
791  {
792  ret[i] = sqrt2 * jacA[i] * jacDerivB[i] * jacC[i] *
793  pow(1.0 - m_eta[2][i], I);
794  }
795  }
796  else
797  {
798  for (int i = 0; i < m_numPoints; ++i)
799  {
800  ret[i] = 2.0 * sqrt2 * jacDerivA[i] * jacB[i] * jacC[i];
801 
802  if (I > 0)
803  {
804  ret[i] *= pow(1.0 - m_eta[2][i], I - 1);
805  }
806  }
807 
808  if (dir == 0)
809  {
810  return ret;
811  }
812 
813  for (int i = 0; i < m_numPoints; ++i)
814  {
815  ret[i] *= 0.5 * (1.0 + m_eta[0][i]);
816 
817  NekDouble tmp = jacDerivC[i] * pow(1.0 - m_eta[2][i], I);
818 
819  if (I > 0)
820  {
821  tmp -= jacC[i] * I * pow(1.0 - m_eta[2][i], I - 1);
822  }
823 
824  ret[i] += sqrt2 * jacA[i] * jacB[i] * tmp;
825  }
826  }
827 
828  return ret;
829 }
830 
831 /**
832  * @brief Construct the nodal utility class for a quadrilateral.
833  *
834  * The constructor of this class sets up the #m_ordering member variable used in
835  * the evaluation of the orthogonal basis.
836  *
837  * @param degree Polynomial order of this nodal quad.
838  * @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
839  * element.
840  * @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
841  * element.
842  */
846  : NodalUtil(degree, 2)
847 {
848  // Set up parent variables.
849  m_numPoints = r.num_elements();
850  m_xi[0] = r;
851  m_xi[1] = s;
852 
853  // Construct a mapping (i,j) -> m from the tensor product space (i,j) to a
854  // single ordering m.
855  for (int j = 0; j <= m_degree; ++j)
856  {
857  for (int i = 0; i <= m_degree; ++i)
858  {
859  m_ordering.push_back(std::make_pair(i,j));
860  }
861  }
862 }
863 
864 /**
865  * @brief Return the value of the modal functions for the quad element at
866  * the nodal points #m_xi for a given mode.
867  *
868  * In a quad, we use the orthogonal basis
869  *
870  * \f[
871  * \psi_{m(ij)} = P^{(0,0)}_i(\xi_1) P_j^{(0,0)}(\xi_2)
872  * \f]
873  *
874  *
875  * @param mode The mode of the orthogonal basis to evaluate.
876  */
878 {
879  std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
880  std::pair<int, int> modes = m_ordering[mode];
881 
882  // Calculate Jacobi polynomials
884  m_numPoints, &m_xi[0][0], &jacobi_i[0], NULL, modes.first, 0.0, 0.0);
886  m_numPoints, &m_xi[1][0], &jacobi_j[0], NULL, modes.second, 0.0, 0.0);
887 
889 
890  for (int i = 0; i < m_numPoints; ++i)
891  {
892  ret[i] = jacobi_i[i] * jacobi_j[i];
893  }
894 
895  return ret;
896 }
897 
898 /**
899  * @brief Return the value of the derivative of the modal functions for the
900  * quadrilateral element at the nodal points #m_xi for a given mode.
901  *
902  * @param mode The mode of the orthogonal basis to evaluate.
903  */
905  const int dir, const int mode)
906 {
907  std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
908  std::vector<NekDouble> jacobi_di(m_numPoints), jacobi_dj(m_numPoints);
909  std::pair<int, int> modes = m_ordering[mode];
910 
911  // Calculate Jacobi polynomials and their derivatives. Note that we use both
912  // jacobfd and jacobd since jacobfd is only valid for derivatives in the
913  // open interval (-1,1).
915  m_numPoints, &m_xi[0][0], &jacobi_i[0], NULL, modes.first, 0.0, 0.0);
917  m_numPoints, &m_xi[1][0], &jacobi_j[0], NULL, modes.second, 0.0, 0.0);
919  m_numPoints, &m_xi[0][0], &jacobi_di[0], modes.first, 0.0, 0.0);
921  m_numPoints, &m_xi[1][0], &jacobi_dj[0], modes.second, 0.0, 0.0);
922 
924 
925  if (dir == 0)
926  {
927  for (int i = 0; i < m_numPoints; ++i)
928  {
929  ret[i] = jacobi_di[i] * jacobi_j[i];
930  }
931  }
932  else
933  {
934  for (int i = 0; i < m_numPoints; ++i)
935  {
936  ret[i] = jacobi_i[i] * jacobi_dj[i];
937  }
938  }
939 
940  return ret;
941 }
942 
943 
944 /**
945  * @brief Construct the nodal utility class for a hexahedron.
946  *
947  * The constructor of this class sets up the #m_ordering member variable used in
948  * the evaluation of the orthogonal basis.
949  *
950  * @param degree Polynomial order of this nodal hexahedron.
951  * @param r \f$ \xi_1 \f$-coordinates of nodal points in the standard
952  * element.
953  * @param s \f$ \xi_2 \f$-coordinates of nodal points in the standard
954  * element.
955  */
960  : NodalUtil(degree, 3)
961 {
962  // Set up parent variables.
963  m_numPoints = r.num_elements();
964  m_xi[0] = r;
965  m_xi[1] = s;
966  m_xi[2] = t;
967 
968  // Construct a mapping (i,j,k) -> m from the tensor product space (i,j,k) to
969  // a single ordering m.
970  for (int k = 0; k <= m_degree; ++k)
971  {
972  for (int j = 0; j <= m_degree; ++j)
973  {
974  for (int i = 0; i <= m_degree; ++i)
975  {
976  m_ordering.push_back(Mode(i, j, k));
977  }
978  }
979  }
980 }
981 
982 /**
983  * @brief Return the value of the modal functions for the hex element at
984  * the nodal points #m_xi for a given mode.
985  *
986  * In a quad, we use the orthogonal basis
987  *
988  * \f[
989  * \psi_{m(ijk)} = P^{(0,0)}_i(\xi_1) P_j^{(0,0)}(\xi_2) P_k^{(0,0)}(\xi_3)
990  * \f]
991  *
992  * @param mode The mode of the orthogonal basis to evaluate.
993  */
995 {
996  std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
997  std::vector<NekDouble> jacobi_k(m_numPoints);
998  Mode modes = m_ordering[mode];
999 
1000  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
1001 
1002  // Calculate Jacobi polynomials
1004  m_numPoints, &m_xi[0][0], &jacobi_i[0], NULL, I, 0.0, 0.0);
1006  m_numPoints, &m_xi[1][0], &jacobi_j[0], NULL, J, 0.0, 0.0);
1008  m_numPoints, &m_xi[2][0], &jacobi_k[0], NULL, K, 0.0, 0.0);
1009 
1011 
1012  for (int i = 0; i < m_numPoints; ++i)
1013  {
1014  ret[i] = jacobi_i[i] * jacobi_j[i] * jacobi_k[i];
1015  }
1016 
1017  return ret;
1018 }
1019 
1021  const int dir, const int mode)
1022 {
1023  std::vector<NekDouble> jacobi_i(m_numPoints), jacobi_j(m_numPoints);
1024  std::vector<NekDouble> jacobi_k(m_numPoints);
1025  std::vector<NekDouble> jacobi_di(m_numPoints), jacobi_dj(m_numPoints);
1026  std::vector<NekDouble> jacobi_dk(m_numPoints);
1027  Mode modes = m_ordering[mode];
1028 
1029  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
1030 
1031  // Calculate Jacobi polynomials and their derivatives. Note that we use both
1032  // jacobfd and jacobd since jacobfd is only valid for derivatives in the
1033  // open interval (-1,1).
1035  m_numPoints, &m_xi[0][0], &jacobi_i[0], NULL, I, 0.0, 0.0);
1037  m_numPoints, &m_xi[1][0], &jacobi_j[0], NULL, J, 0.0, 0.0);
1039  m_numPoints, &m_xi[2][0], &jacobi_k[0], NULL, K, 0.0, 0.0);
1041  m_numPoints, &m_xi[0][0], &jacobi_di[0], I, 0.0, 0.0);
1043  m_numPoints, &m_xi[1][0], &jacobi_dj[0], J, 0.0, 0.0);
1045  m_numPoints, &m_xi[2][0], &jacobi_dk[0], K, 0.0, 0.0);
1046 
1048 
1049  if (dir == 0)
1050  {
1051  for (int i = 0; i < m_numPoints; ++i)
1052  {
1053  ret[i] = jacobi_di[i] * jacobi_j[i] * jacobi_k[i];
1054  }
1055  }
1056  else if (dir == 1)
1057  {
1058  for (int i = 0; i < m_numPoints; ++i)
1059  {
1060  ret[i] = jacobi_dj[i] * jacobi_i[i] * jacobi_k[i];
1061  }
1062  }
1063  else
1064  {
1065  for (int i = 0; i < m_numPoints; ++i)
1066  {
1067  ret[i] = jacobi_i[i] * jacobi_j[i] * jacobi_dk[i];
1068  }
1069  }
1070 
1071  return ret;
1072 }
1073 
1074 }
1075 }
1076 
std::vector< Mode > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:232
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)
Return the value of the derivative of the modal functions for the tetrahedral element at the nodal po...
Definition: NodalUtil.cpp:537
A class to assist in the construction of nodal simplex and hybrid elements in two and three dimension...
Definition: NodalUtil.h:84
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)
Return the value of the modal functions for the quad element at the nodal points m_xi for a given mod...
Definition: NodalUtil.cpp:877
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
boost::tuple< int, int, int > Mode
Definition: NodalUtil.h:355
NodalUtilPrism(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s, Array< OneD, NekDouble > t)
Construct the nodal utility class for a prism.
Definition: NodalUtil.cpp:655
Array< OneD, Array< OneD, NekDouble > > m_eta
Collapsed coordinates of the nodal points.
Definition: NodalUtil.h:285
Array< OneD, Array< OneD, NekDouble > > m_eta
Collapsed coordinates of the nodal points.
Definition: NodalUtil.h:236
SharedMatrix GetVandermondeForDeriv(int dir)
Return the Vandermonde matrix of the derivative of the basis functions for the nodal distribution...
Definition: NodalUtil.cpp:135
SharedMatrix GetInterpolationMatrix(Array< OneD, Array< OneD, NekDouble > > &xi)
Construct the interpolation matrix used to evaluate the basis at the points xi inside the element...
Definition: NodalUtil.cpp:202
boost::tuple< int, int, int > Mode
Definition: NodalUtil.h:217
std::vector< std::pair< int, int > > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:184
int m_degree
Degree of the nodal element.
Definition: NodalUtil.h:108
boost::shared_ptr< NekMatrix< NekDouble > > SharedMatrix
Definition: NodalUtil.h:56
static const NekDouble kNekZeroTol
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:2151
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)
Return the value of the modal functions for the triangular element at the nodal points m_xi for a giv...
Definition: NodalUtil.cpp:295
Array< OneD, Array< OneD, NekDouble > > m_eta
Collapsed coordinates of the nodal points.
Definition: NodalUtil.h:187
SharedMatrix GetDerivMatrix(int dir)
Return the derivative matrix for the nodal distribution.
Definition: NodalUtil.cpp:168
std::vector< Mode > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:370
double NekDouble
NekVector< NekDouble > GetWeights()
Obtain the integration weights for the given nodal distribution.
Definition: NodalUtil.cpp:66
virtual boost::shared_ptr< NodalUtil > v_CreateUtil(Array< OneD, Array< OneD, NekDouble > > &xi)=0
Construct a NodalUtil object of the appropriate element type for a given set of points.
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)
Return the value of the derivative of the modal functions for the quadrilateral element at the nodal ...
Definition: NodalUtil.cpp:904
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)
Return the value of the modal functions for the prismatic element at the nodal points m_xi for a give...
Definition: NodalUtil.cpp:718
SharedMatrix GetVandermonde()
Return the Vandermonde matrix for the nodal distribution.
Definition: NodalUtil.cpp:102
NodalUtilHex(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s, Array< OneD, NekDouble > t)
Construct the nodal utility class for a hexahedron.
Definition: NodalUtil.cpp:956
NodalUtilQuad(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s)
Construct the nodal utility class for a quadrilateral.
Definition: NodalUtil.cpp:843
Array< OneD, Array< OneD, NekDouble > > m_xi
Coordinates of the nodal points defining the basis.
Definition: NodalUtil.h:112
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)
Return the values of the derivative of the orthogonal basis at the nodal points for a given mode...
Definition: NodalUtil.cpp:1020
std::vector< Mode > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:281
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)=0
Return the values of the orthogonal basis at the nodal points for a given mode.
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)
Return the value of the derivative of the modal functions for the prismatic element at the nodal poin...
Definition: NodalUtil.cpp:760
int m_numPoints
Total number of nodal points.
Definition: NodalUtil.h:110
NodalUtilTriangle(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s)
Construct the nodal utility class for a triangle.
Definition: NodalUtil.cpp:239
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)
Return the value of the modal functions for the hex element at the nodal points m_xi for a given mode...
Definition: NodalUtil.cpp:994
boost::tuple< int, int, int > Mode
Definition: NodalUtil.h:266
std::vector< std::pair< int, int > > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:326
virtual int v_NumModes()=0
Calculate the number of degrees of freedom for this element.
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)
Return the value of the derivative of the modal functions for the triangular element at the nodal poi...
Definition: NodalUtil.cpp:333
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)=0
Return the values of the derivative of the orthogonal basis at the nodal points for a given mode...
RawType< VectorType >::type SolveTranspose(const VectorType &b)
Definition: NekLinSys.hpp:421
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)
Return the value of the modal functions for the tetrahedral element at the nodal points m_xi for a gi...
Definition: NodalUtil.cpp:495
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1951
virtual NekDouble v_ModeZeroIntegral()=0
Return the value of the integral of the zero-th mode for this element.
NodalUtilTetrahedron(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s, Array< OneD, NekDouble > t)
Construct the nodal utility class for a tetrahedron.
Definition: NodalUtil.cpp:420