Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdExpansion.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Stdexpansion.h
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: Class definition StdExpansion which is the base class
33 // to all expansion shapes
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #ifndef NEKTAR_LIB_STDREGIONS_STANDARDEXPANSION_H
38 #define NEKTAR_LIB_STDREGIONS_STANDARDEXPANSION_H
39 
40 #include <fstream>
41 #include <vector>
42 
47 #include <StdRegions/IndexMapKey.h>
49 #include <boost/enable_shared_from_this.hpp>
50 namespace Nektar { namespace LocalRegions { class MatrixKey; class Expansion; } }
51 
52 
53 namespace Nektar
54 {
55  namespace StdRegions
56  {
57 
58  class StdExpansion1D;
60 
61  typedef Array<OneD, Array<OneD, NekDouble> > NormalVector;
62 
63  /** \brief The base class for all shapes
64  *
65  * This is the lowest level basic class for all shapes and so
66  * contains the definition of common data and common routine to all
67  * elements
68  */
69  class StdExpansion : public boost::enable_shared_from_this<StdExpansion>
70  {
71  public:
72 
73  /** \brief Default Constructor */
75 
76  /** \brief Constructor */
77  STD_REGIONS_EXPORT StdExpansion(const int numcoeffs, const int numbases,
81 
82 
83  /** \brief Copy Constructor */
84  STD_REGIONS_EXPORT StdExpansion(const StdExpansion &T);
85 
86  /** \brief Destructor */
88 
89 
90  // Standard Expansion Routines Applicable Regardless of Region
91 
92  /** \brief This function returns the number of 1D bases used in
93  * the expansion
94  *
95  * \return returns the number of 1D bases used in the expansion,
96  * which is equal to number dimension of the expansion
97  */
98  inline int GetNumBases() const
99  {
100  return m_base.num_elements();
101  }
102 
103  /** \brief This function gets the shared point to basis
104  *
105  * \return returns the shared pointer to the bases
106  */
107  inline const Array<OneD, const LibUtilities::BasisSharedPtr>& GetBase() const
108  {
109  return(m_base);
110  }
111 
112  /** \brief This function gets the shared point to basis in
113  * the \a dir direction
114  *
115  * \return returns the shared pointer to the basis in
116  * directin \a dir
117  */
118  inline const LibUtilities::BasisSharedPtr& GetBasis(int dir) const
119  {
120  ASSERTL1(dir < m_base.num_elements(),
121  "dir is larger than number of bases");
122  return(m_base[dir]);
123  }
124 
125  /** \brief This function returns the total number of coefficients
126  * used in the expansion
127  *
128  * \return returns the total number of coefficients (which is
129  * equivalent to the total number of modes) used in the expansion
130  */
131  inline int GetNcoeffs(void) const
132  {
133  return(m_ncoeffs);
134  }
135 
136  /** \brief This function returns the total number of quadrature
137  * points used in the element
138  *
139  * \return returns the total number of quadrature points
140  */
141  inline int GetTotPoints() const
142  {
143  int i;
144  int nqtot = 1;
145 
146  for(i=0; i < m_base.num_elements(); ++i)
147  {
148  nqtot *= m_base[i]->GetNumPoints();
149  }
150 
151  return nqtot;
152  }
153 
154 
155  /** \brief This function returns the type of basis used in the \a dir
156  * direction
157  *
158  * The different types of bases implemented in the code are defined
159  * in the LibUtilities::BasisType enumeration list. As a result, the
160  * function will return one of the types of this enumeration list.
161  *
162  * \param dir the direction
163  * \return returns the type of basis used in the \a dir direction
164  */
165  inline LibUtilities::BasisType GetBasisType(const int dir) const
166  {
167  ASSERTL1(dir < m_base.num_elements(), "dir is larger than m_numbases");
168  return(m_base[dir]->GetBasisType());
169  }
170 
171  /** \brief This function returns the number of expansion modes
172  * in the \a dir direction
173  *
174  * \param dir the direction
175  * \return returns the number of expansion modes in the \a dir
176  * direction
177  */
178  inline int GetBasisNumModes(const int dir) const
179  {
180  ASSERTL1(dir < m_base.num_elements(),"dir is larger than m_numbases");
181  return(m_base[dir]->GetNumModes());
182  }
183 
184 
185  /** \brief This function returns the maximum number of
186  * expansion modes over all local directions
187  *
188  * \return returns the maximum number of expansion modes
189  * over all local directions
190  */
191  inline int EvalBasisNumModesMax(void) const
192  {
193  int i;
194  int returnval = 0;
195 
196  for(i = 0; i < m_base.num_elements(); ++i)
197  {
198  returnval = max(returnval, m_base[i]->GetNumModes());
199  }
200 
201  return returnval;
202  }
203 
204  /** \brief This function returns the type of quadrature points used
205  * in the \a dir direction
206  *
207  * The different types of quadrature points implemented in the code
208  * are defined in the LibUtilities::PointsType enumeration list.
209  * As a result, the function will return one of the types of this
210  * enumeration list.
211  *
212  * \param dir the direction
213  * \return returns the type of quadrature points used in the \a dir
214  * direction
215  */
216  inline LibUtilities::PointsType GetPointsType(const int dir) const
217  {
218  ASSERTL1(dir < m_base.num_elements(), "dir is larger than m_numbases");
219  return(m_base[dir]->GetPointsType());
220  }
221 
222  /** \brief This function returns the number of quadrature points
223  * in the \a dir direction
224  *
225  * \param dir the direction
226  * \return returns the number of quadrature points in the \a dir
227  * direction
228  */
229  inline int GetNumPoints(const int dir) const
230  {
231  ASSERTL1(dir < m_base.num_elements() || dir == 0,
232  "dir is larger than m_numbases");
233  return(m_base.num_elements() > 0 ? m_base[dir]->GetNumPoints() : 1);
234  }
235 
236  /** \brief This function returns a pointer to the array containing
237  * the quadrature points in \a dir direction
238  *
239  * \param dir the direction
240  * \return returns a pointer to the array containing
241  * the quadrature points in \a dir direction
242  */
243  inline const Array<OneD, const NekDouble>& GetPoints(const int dir) const
244  {
245  return m_base[dir]->GetZ();
246  }
247 
248 
249  // Wrappers around virtual Functions
250 
251  /** \brief This function returns the number of vertices of the
252  * expansion domain
253  *
254  * This function is a wrapper around the virtual function
255  * \a v_GetNverts()
256  *
257  * \return returns the number of vertices of the expansion domain
258  */
259  int GetNverts() const
260  {
261  return v_GetNverts();
262  }
263 
264  /** \brief This function returns the number of edges of the
265  * expansion domain
266  *
267  * This function is a wrapper around the virtual function
268  * \a v_GetNedges()
269  *
270  * \return returns the number of edges of the expansion domain
271  */
272  int GetNedges() const
273  {
274  return v_GetNedges();
275  }
276 
277  /** \brief This function returns the number of expansion coefficients
278  * belonging to the \a i-th edge
279  *
280  * This function is a wrapper around the virtual function
281  * \a v_GetEdgeNcoeffs()
282  *
283  * \param i specifies which edge
284  * \return returns the number of expansion coefficients belonging to
285  * the \a i-th edge
286  */
287  int GetEdgeNcoeffs(const int i) const
288  {
289  return v_GetEdgeNcoeffs(i);
290  }
291 
292 
294  {
295  return v_GetTotalEdgeIntNcoeffs();
296  }
297 
298  /** \brief This function returns the number of quadrature points
299  * belonging to the \a i-th edge
300  *
301  * This function is a wrapper around the virtual function
302  * \a v_GetEdgeNumPoints()
303  *
304  * \param i specifies which edge
305  * \return returns the number of expansion coefficients belonging to
306  * the \a i-th edge
307  */
308  int GetEdgeNumPoints(const int i) const
309  {
310  return v_GetEdgeNumPoints(i);
311  }
312 
313 
314  int DetCartesianDirOfEdge(const int edge)
315  {
316  return v_DetCartesianDirOfEdge(edge);
317  }
318 
319  const LibUtilities::BasisKey DetEdgeBasisKey(const int i) const
320  {
321  return v_DetEdgeBasisKey(i);
322  }
323 
324  const LibUtilities::BasisKey DetFaceBasisKey(const int i, const int k) const
325  {
326  return v_DetFaceBasisKey(i, k);
327  }
328  /**
329  * \brief This function returns the number of quadrature points
330  * belonging to the \a i-th face.
331  *
332  * This function is a wrapper around the virtual function \a
333  * v_GetFaceNcoeffs()
334  *
335  * \param i specifies which face
336  * \return returns the number of expansion coefficients belonging to
337  * the \a i-th face
338  */
339  int GetFaceNumPoints(const int i) const
340  {
341  return v_GetFaceNumPoints(i);
342  }
343 
344  /** \brief This function returns the number of expansion coefficients
345  * belonging to the \a i-th face
346  *
347  * This function is a wrapper around the virtual function
348  * \a v_GetFaceNcoeffs()
349  *
350  * \param i specifies which face
351  * \return returns the number of expansion coefficients belonging to
352  * the \a i-th face
353  */
354  int GetFaceNcoeffs(const int i) const
355  {
356  return v_GetFaceNcoeffs(i);
357  }
358 
359  int GetFaceIntNcoeffs(const int i) const
360  {
361  return v_GetFaceIntNcoeffs(i);
362  }
363 
365  {
366  return v_GetTotalFaceIntNcoeffs();
367  }
368 
369  LibUtilities::PointsKey GetFacePointsKey(const int i, const int j) const
370  {
371  return v_GetFacePointsKey(i, j);
372  }
373 
374  int NumBndryCoeffs(void) const
375  {
376  return v_NumBndryCoeffs();
377  }
378 
379  int NumDGBndryCoeffs(void) const
380  {
381  return v_NumDGBndryCoeffs();
382  }
383 
384  /** \brief This function returns the type of expansion basis on the
385  * \a i-th edge
386  *
387  * This function is a wrapper around the virtual function
388  * \a v_GetEdgeBasisType()
389  *
390  * The different types of bases implemented in the code are defined
391  * in the LibUtilities::BasisType enumeration list. As a result, the
392  * function will return one of the types of this enumeration list.
393  *
394  * \param i specifies which edge
395  * \return returns the expansion basis on the \a i-th edge
396  */
398  {
399  return v_GetEdgeBasisType(i);
400  }
401 
402 
403  /** \brief This function returns the number of faces of the
404  * expansion domain
405  *
406  * This function is a wrapper around the virtual function
407  * \a v_GetNFaces()
408  *
409  * \return returns the number of faces of the expansion domain
410  */
411  int GetNfaces() const
412  {
413  return v_GetNfaces();
414  }
415 
416  /** \brief This function returns the shape of the expansion domain
417  *
418  * This function is a wrapper around the virtual function
419  * \a v_DetShapeType()
420  *
421  * The different shape types implemented in the code are defined
422  * in the ::ShapeType enumeration list. As a result, the
423  * function will return one of the types of this enumeration list.
424  *
425  * \return returns the shape of the expansion domain
426  */
428  {
429  return v_DetShapeType();
430  }
431 
432  int GetShapeDimension() const
433  {
434  return v_GetShapeDimension();
435  }
436 
438  {
440  }
441 
442 
443  /** \brief This function performs the Backward transformation from
444  * coefficient space to physical space
445  *
446  * This function is a wrapper around the virtual function
447  * \a v_BwdTrans()
448  *
449  * Based on the expansion coefficients, this function evaluates the
450  * expansion at the quadrature points. This is equivalent to the
451  * operation \f[ u(\xi_{1i}) =
452  * \sum_{p=0}^{P-1} \hat{u}_p \phi_p(\xi_{1i}) \f] which can be
453  * evaluated as \f$ {\bf u} = {\bf B}^T {\bf \hat{u}} \f$ with
454  * \f${\bf B}[i][j] = \phi_i(\xi_{j})\f$
455  *
456  * This function requires that the coefficient array
457  * \f$\mathbf{\hat{u}}\f$ provided as \a inarray.
458  *
459  * The resulting array
460  * \f$\mathbf{u}[m]=u(\mathbf{\xi}_m)\f$ containing the
461  * expansion evaluated at the quadrature points, is stored
462  * in the \a outarray.
463  *
464  * \param inarray contains the values of the expansion
465  * coefficients (input of the function)
466  *
467  * \param outarray contains the values of the expansion evaluated
468  * at the quadrature points (output of the function)
469  */
470 
471  void BwdTrans (const Array<OneD, const NekDouble>& inarray,
472  Array<OneD, NekDouble> &outarray)
473  {
474  v_BwdTrans (inarray, outarray);
475  }
476 
477  /**
478  * @brief This function performs the Forward transformation from
479  * physical space to coefficient space.
480  */
481  inline void FwdTrans (const Array<OneD, const NekDouble>& inarray,
482  Array<OneD, NekDouble> &outarray);
483 
484  void FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
485  Array<OneD, NekDouble> &outarray)
486  {
487  v_FwdTrans_BndConstrained(inarray,outarray);
488  }
489 
490  /** \brief This function integrates the specified function over the
491  * domain
492  *
493  * This function is a wrapper around the virtual function
494  * \a v_Integral()
495  *
496  * Based on the values of the function evaluated at the quadrature
497  * points (which are stored in \a inarray), this function calculates
498  * the integral of this function over the domain. This is
499  * equivalent to the numerical evaluation of the operation
500  * \f[ I=\int u(\mathbf{\xi})d \mathbf{\xi}\f]
501  *
502  * \param inarray values of the function to be integrated evaluated
503  * at the quadrature points (i.e.
504  * \a inarray[m]=\f$u(\mathbf{\xi}_m)\f$)
505  * \return returns the value of the calculated integral
506  *
507  * Inputs:\n
508 
509  - \a inarray: definition of function to be returned at quadrature point
510  of expansion.
511 
512  Outputs:\n
513 
514  - returns \f$\int^1_{-1}\int^1_{-1} u(\xi_1, \xi_2) J[i,j] d
515  \xi_1 d \xi_2 \f$ where \f$inarray[i,j] =
516  u(\xi_{1i},\xi_{2j}) \f$ and \f$ J[i,j] \f$ is the
517  Jacobian evaluated at the quadrature point.
518  *
519  */
520  NekDouble Integral(const Array<OneD, const NekDouble>& inarray )
521  {
522  return v_Integral(inarray);
523  }
524 
525  /** \brief This function fills the array \a outarray with the
526  * \a mode-th mode of the expansion
527  *
528  * This function is a wrapper around the virtual function
529  * \a v_FillMode()
530  *
531  * The requested mode is evaluated at the quadrature points
532  *
533  * \param mode the mode that should be filled
534  * \param outarray contains the values of the \a mode-th mode of the
535  * expansion evaluated at the quadrature points (output of the
536  * function)
537  */
538  void FillMode(const int mode, Array<OneD, NekDouble> &outarray)
539  {
540  v_FillMode(mode, outarray);
541  }
542 
543  /** \brief this function calculates the inner product of a given
544  * function \a f with the different modes of the expansion
545  *
546  * This function is a wrapper around the virtual function
547  * \a v_IProductWRTBase()
548  *
549  * This is equivalent to the numerical evaluation of
550  * \f[ I[p] = \int \phi_p(\mathbf{x}) f(\mathbf{x}) d\mathbf{x}\f]
551  * \f$ \begin{array}{rcl} I_{pq} = (\phi_q \phi_q, u) & = &
552  \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \phi_p(\xi_{0,i})
553  \phi_q(\xi_{1,j}) w^0_i w^1_j u(\xi_{0,i} \xi_{1,j})
554  J_{i,j}\\ & = & \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i})
555  \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j}
556  J_{i,j} \end{array} \f$
557 
558  where
559 
560  \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
561 
562  which can be implemented as
563 
564  \f$ f_{qi} = \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j} =
565  {\bf B_1 U} \f$
566  \f$ I_{pq} = \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i}) f_{qi} =
567  {\bf B_0 F} \f$
568  *
569  * \param inarray contains the values of the function \a f
570  * evaluated at the quadrature points
571  * \param outarray contains the values of the inner product of \a f
572  * with the different modes, i.e. \f$ outarray[p] = I[p]\f$
573  * (output of the function)
574  */
575  void IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
576  Array<OneD, NekDouble> &outarray)
577  {
578  v_IProductWRTBase(inarray, outarray);
579  }
580 
582  const Array<OneD, const NekDouble>& base,
583  const Array<OneD, const NekDouble>& inarray,
584  Array<OneD, NekDouble> &outarray,
585  int coll_check)
586  {
587  v_IProductWRTBase(base, inarray, outarray, coll_check);
588  }
589 
590 
591  void IProductWRTDerivBase(const int dir,
592  const Array<OneD, const NekDouble>& inarray,
593  Array<OneD, NekDouble> &outarray)
594  {
595  v_IProductWRTDerivBase(dir,inarray, outarray);
596  }
597 
598  /// \brief Get the element id of this expansion when used
599  /// in a list by returning value of #m_elmt_id
600  inline int GetElmtId()
601  {
602  return m_elmt_id;
603  }
604 
605 
606  /// \brief Set the element id of this expansion when used
607  /// in a list by returning value of #m_elmt_id
608  inline void SetElmtId(const int id)
609  {
610  m_elmt_id = id;
611  }
612 
613  /** \brief this function returns the physical coordinates of the
614  * quadrature points of the expansion
615  *
616  * This function is a wrapper around the virtual function
617  * \a v_GetCoords()
618  *
619  * \param coords an array containing the coordinates of the
620  * quadrature points (output of the function)
621  */
622  void GetCoords(Array<OneD, NekDouble> &coords_1,
623  Array<OneD, NekDouble> &coords_2 = NullNekDouble1DArray,
624  Array<OneD, NekDouble> &coords_3 = NullNekDouble1DArray)
625  {
626  v_GetCoords(coords_1,coords_2,coords_3);
627  }
628 
629  /** \brief given the coordinates of a point of the element in the
630  * local collapsed coordinate system, this function calculates the
631  * physical coordinates of the point
632  *
633  * This function is a wrapper around the virtual function
634  * \a v_GetCoord()
635  *
636  * \param Lcoords the coordinates in the local collapsed
637  * coordinate system
638  * \param coords the physical coordinates (output of the function)
639  */
640  void GetCoord(const Array<OneD, const NekDouble>& Lcoord,
641  Array<OneD, NekDouble> &coord)
642  {
643  v_GetCoord(Lcoord, coord);
644  }
645 
647  {
648  return m_stdMatrixManager[mkey];
649  }
650 
652  {
653  return m_stdStaticCondMatrixManager[mkey];
654  }
655 
657  {
658  return m_IndexMapManager[ikey];
659  }
660 
661  const Array<OneD, const NekDouble>& GetPhysNormals(void)
662  {
663  return v_GetPhysNormals();
664  }
665 
666  void SetPhysNormals(Array<OneD, const NekDouble> &normal)
667  {
668  v_SetPhysNormals(normal);
669  }
670 
671  STD_REGIONS_EXPORT virtual void SetUpPhysNormals(const int edge);
672 
673  void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, Array< OneD, NekDouble> &outarray)
674  {
675  v_NormVectorIProductWRTBase(Fx,Fy,outarray);
676  }
677 
678  void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, const Array<OneD, const NekDouble> &Fz, Array< OneD, NekDouble> &outarray)
679  {
680  v_NormVectorIProductWRTBase(Fx,Fy,Fz,outarray);
681  }
682 
684  {
685  return v_GetLocStaticCondMatrix(mkey);
686  }
687 
689  {
690  return v_DropLocStaticCondMatrix(mkey);
691  }
692 
694  {
695  return v_GetFaceOrient(face);
696  }
697 
699  {
700  return v_GetEorient(edge);
701  }
702 
704  {
705  return v_GetPorient(point);
706  }
707 
709  {
710  return v_GetCartesianEorient(edge);
711  }
712 
714  Array<OneD, NekDouble> &coeffs,
716  {
717  v_SetCoeffsToOrientation(coeffs, dir);
718  }
719 
722  Array<OneD, const NekDouble> &inarray,
723  Array<OneD, NekDouble> &outarray)
724  {
725  v_SetCoeffsToOrientation(dir,inarray,outarray);
726  }
727 
728  int CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset)
729  {
730  return v_CalcNumberOfCoefficients(nummodes,modes_offset);
731  }
732 
733  void ExtractDataToCoeffs(const NekDouble *data,
734  const std::vector<unsigned int > &nummodes,
735  const int nmodes_offset,
736  NekDouble *coeffs)
737  {
738  v_ExtractDataToCoeffs(data,nummodes,nmodes_offset,coeffs);
739  }
740 
741  // virtual functions related to LocalRegions
743  const Array<OneD, const NekDouble> &Lcoord,
744  const Array<OneD, const NekDouble> &physvals);
745 
746 
747  STD_REGIONS_EXPORT void AddEdgeNormBoundaryInt(const int edge,
748  boost::shared_ptr<StdExpansion> &EdgeExp,
749  const Array<OneD, const NekDouble> &Fx,
750  const Array<OneD, const NekDouble> &Fy,
751  Array<OneD, NekDouble> &outarray);
752 
753  STD_REGIONS_EXPORT void AddEdgeNormBoundaryInt(const int edge,
754  boost::shared_ptr<StdExpansion> &EdgeExp,
755  const Array<OneD, const NekDouble> &Fn,
756  Array<OneD, NekDouble> &outarray);
757 
758  STD_REGIONS_EXPORT void AddEdgeNormBoundaryBiInt(const int edge,
759  boost::shared_ptr<StdExpansion> &EdgeExp,
760  const Array<OneD, const NekDouble> &Fwd,
761  const Array<OneD, const NekDouble> &Bwd,
762  Array<OneD, NekDouble> &outarray);
763 
764  STD_REGIONS_EXPORT void AddFaceNormBoundaryInt(const int face,
765  boost::shared_ptr<StdExpansion> &FaceExp,
766  const Array<OneD, const NekDouble> &Fn,
767  Array<OneD, NekDouble> &outarray);
768 
770  {
771  return v_GetCoordim();
772  }
773 
774  void GetBoundaryMap(Array<OneD, unsigned int> &outarray)
775  {
776  v_GetBoundaryMap(outarray);
777  }
778 
779  void GetInteriorMap(Array<OneD, unsigned int> &outarray)
780  {
781  v_GetInteriorMap(outarray);
782  }
783 
784  int GetVertexMap(const int localVertexId,
785  bool useCoeffPacking = false)
786  {
787  return v_GetVertexMap(localVertexId,useCoeffPacking);
788  }
789 
790  void GetEdgeInteriorMap(const int eid, const Orientation edgeOrient,
791  Array<OneD, unsigned int> &maparray,
792  Array<OneD, int> &signarray)
793  {
794  v_GetEdgeInteriorMap(eid,edgeOrient,maparray,signarray);
795  }
796 
797  void GetFaceInteriorMap(const int fid, const Orientation faceOrient,
798  Array<OneD, unsigned int> &maparray,
799  Array<OneD, int> &signarray)
800  {
801  v_GetFaceInteriorMap(fid,faceOrient,maparray,signarray);
802  }
803 
804  void GetEdgeToElementMap(const int eid, const Orientation edgeOrient,
805  Array<OneD, unsigned int> &maparray,
806  Array<OneD, int> &signarray)
807  {
808  v_GetEdgeToElementMap(eid,edgeOrient,maparray,signarray);
809  }
810 
811  void GetFaceToElementMap(const int fid, const Orientation faceOrient,
812  Array<OneD, unsigned int> &maparray,
813  Array<OneD, int> &signarray,
814  int nummodesA = -1, int nummodesB = -1)
815  {
816  v_GetFaceToElementMap(fid,faceOrient,maparray,signarray,
817  nummodesA,nummodesB);
818  }
819 
820 
821  /**
822  * @brief Extract the physical values along edge \a edge from \a
823  * inarray into \a outarray following the local edge orientation
824  * and point distribution defined by defined in \a EdgeExp.
825  */
826 
827  void GetEdgePhysVals(const int edge, const Array<OneD,
828  const NekDouble> &inarray,
829  Array<OneD,NekDouble> &outarray)
830  {
831  v_GetEdgePhysVals(edge,inarray,outarray);
832  }
833 
834  void GetEdgePhysVals(const int edge,
835  const boost::shared_ptr<StdExpansion> &EdgeExp,
836  const Array<OneD, const NekDouble> &inarray,
837  Array<OneD,NekDouble> &outarray)
838  {
839  v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
840  }
841 
842  void GetTracePhysVals(const int edge, const boost::shared_ptr<StdExpansion> &EdgeExp, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
843  {
844  v_GetTracePhysVals(edge,EdgeExp,inarray,outarray);
845  }
846 
847  void GetVertexPhysVals(const int vertex,
848  const Array<OneD, const NekDouble> &inarray,
849  NekDouble &outarray)
850  {
851  v_GetVertexPhysVals(vertex, inarray, outarray);
852  }
853 
854  void GetEdgeInterpVals(const int edge,const Array<OneD,
855  const NekDouble> &inarray,
856  Array<OneD,NekDouble> &outarray)
857  {
858  v_GetEdgeInterpVals(edge, inarray, outarray);
859  }
860 
861  /**
862  * @brief Extract the metric factors to compute the contravariant
863  * fluxes along edge \a edge and stores them into \a outarray
864  * following the local edge orientation (i.e. anticlockwise
865  * convention).
866  */
868  const int edge,
869  Array<OneD, NekDouble> &outarray)
870  {
871  v_GetEdgeQFactors(edge, outarray);
872  }
873 
874 
875 
877  const int face,
878  const boost::shared_ptr<StdExpansion> &FaceExp,
879  const Array<OneD, const NekDouble> &inarray,
880  Array<OneD, NekDouble> &outarray,
882  {
883  v_GetFacePhysVals(face, FaceExp, inarray, outarray, orient);
884  }
885 
887  const Array<OneD, const NekDouble> &inarray,
888  Array<OneD, NekDouble> &outarray)
889  {
890  v_MultiplyByQuadratureMetric(inarray, outarray);
891  }
892 
894  const Array<OneD, const NekDouble> &inarray,
895  Array<OneD, NekDouble> & outarray)
896  {
897  v_MultiplyByStdQuadratureMetric(inarray, outarray);
898  }
899 
900  // Matrix Routines
901 
902  /** \brief this function generates the mass matrix
903  * \f$\mathbf{M}[i][j] =
904  * \int \phi_i(\mathbf{x}) \phi_j(\mathbf{x}) d\mathbf{x}\f$
905  *
906  * \return returns the mass matrix
907  */
908 
910 
911  STD_REGIONS_EXPORT void GeneralMatrixOp(const Array<OneD, const NekDouble> &inarray,
912  Array<OneD,NekDouble> &outarray,
913  const StdMatrixKey &mkey);
914 
915  void MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
916  Array<OneD,NekDouble> &outarray,
917  const StdMatrixKey &mkey)
918  {
919  v_MassMatrixOp(inarray,outarray,mkey);
920  }
921 
922  void LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
923  Array<OneD,NekDouble> &outarray,
924  const StdMatrixKey &mkey)
925  {
926  v_LaplacianMatrixOp(inarray,outarray,mkey);
927  }
928 
929  void ReduceOrderCoeffs(int numMin,
930  const Array<OneD, const NekDouble> &inarray,
931  Array<OneD,NekDouble> &outarray)
932  {
933  v_ReduceOrderCoeffs(numMin,inarray,outarray);
934  }
935 
936  void SVVLaplacianFilter(Array<OneD,NekDouble> &array,
937  const StdMatrixKey &mkey)
938  {
939  v_SVVLaplacianFilter(array,mkey);
940  }
941 
942  void LaplacianMatrixOp(const int k1, const int k2,
943  const Array<OneD, const NekDouble> &inarray,
944  Array<OneD,NekDouble> &outarray,
945  const StdMatrixKey &mkey)
946  {
947  v_LaplacianMatrixOp(k1,k2,inarray,outarray,mkey);
948  }
949 
950  void WeakDerivMatrixOp(const int i,
951  const Array<OneD, const NekDouble> &inarray,
952  Array<OneD,NekDouble> &outarray,
953  const StdMatrixKey &mkey)
954  {
955  v_WeakDerivMatrixOp(i,inarray,outarray,mkey);
956  }
957 
958  void WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
959  Array<OneD,NekDouble> &outarray,
960  const StdMatrixKey &mkey)
961  {
962  v_WeakDirectionalDerivMatrixOp(inarray,outarray,mkey);
963  }
964 
965  void MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
966  Array<OneD,NekDouble> &outarray,
967  const StdMatrixKey &mkey)
968  {
969  v_MassLevelCurvatureMatrixOp(inarray,outarray,mkey);
970  }
971 
972  void LinearAdvectionDiffusionReactionMatrixOp(const Array<OneD, const NekDouble> &inarray,
973  Array<OneD,NekDouble> &outarray,
974  const StdMatrixKey &mkey,
975  bool addDiffusionTerm = true)
976  {
977  v_LinearAdvectionDiffusionReactionMatrixOp(inarray,outarray,mkey,addDiffusionTerm);
978  }
979 
980  /**
981  * @param inarray Input array @f$ \mathbf{u} @f$.
982  * @param outarray Output array @f$ \boldsymbol{\nabla^2u}
983  * + \lambda \boldsymbol{u} @f$.
984  * @param mkey
985  */
986  void HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
987  Array<OneD,NekDouble> &outarray,
988  const StdMatrixKey &mkey)
989  {
990  v_HelmholtzMatrixOp(inarray,outarray,mkey);
991  }
992 
994  {
995  return v_GenMatrix(mkey);
996  }
997 
998  void PhysDeriv (const Array<OneD, const NekDouble>& inarray,
999  Array<OneD, NekDouble> &out_d0,
1000  Array<OneD, NekDouble> &out_d1 = NullNekDouble1DArray,
1001  Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray)
1002  {
1003  v_PhysDeriv (inarray, out_d0, out_d1, out_d2);
1004  }
1005 
1006  void PhysDeriv(const int dir,
1007  const Array<OneD, const NekDouble>& inarray,
1008  Array<OneD, NekDouble> &outarray)
1009  {
1010  v_PhysDeriv (dir, inarray, outarray);
1011  }
1012 
1013  void PhysDeriv_s(const Array<OneD, const NekDouble>& inarray,
1014  Array<OneD, NekDouble> &out_ds)
1015  {
1016  v_PhysDeriv_s(inarray,out_ds);
1017  }
1018 
1019  void PhysDeriv_n(const Array<OneD, const NekDouble>& inarray,
1020  Array<OneD, NekDouble>& out_dn)
1021  {
1022  v_PhysDeriv_n(inarray,out_dn);
1023  }
1024 
1025  void PhysDirectionalDeriv(const Array<OneD, const NekDouble>& inarray,
1026  const Array<OneD, const NekDouble>& direction,
1027  Array<OneD, NekDouble> &outarray)
1028  {
1029  v_PhysDirectionalDeriv (inarray, direction, outarray);
1030  }
1031 
1032  void StdPhysDeriv(const Array<OneD, const NekDouble>& inarray,
1033  Array<OneD, NekDouble> &out_d0,
1034  Array<OneD, NekDouble> &out_d1 = NullNekDouble1DArray,
1035  Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray)
1036  {
1037  v_StdPhysDeriv(inarray, out_d0, out_d1, out_d2);
1038  }
1039 
1040  void StdPhysDeriv (const int dir,
1041  const Array<OneD, const NekDouble>& inarray,
1042  Array<OneD, NekDouble> &outarray)
1043  {
1044  v_StdPhysDeriv(dir,inarray,outarray);
1045  }
1046 
1047  void AddRobinMassMatrix(const int edgeid, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
1048  {
1049  v_AddRobinMassMatrix(edgeid,primCoeffs,inoutmat);
1050  }
1051 
1052  void AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble> &primCoeffs, Array<OneD, NekDouble> &coeffs)
1053  {
1054  v_AddRobinEdgeContribution(edgeid, primCoeffs, coeffs);
1055  }
1056 
1057  void DGDeriv(const int dir,
1058  const Array<OneD, const NekDouble>& inarray,
1059  Array<OneD, boost::shared_ptr< StdExpansion > > &EdgeExp,
1060  Array<OneD, Array<OneD, NekDouble> > &coeffs,
1061  Array<OneD, NekDouble> &outarray)
1062  {
1063  v_DGDeriv (dir, inarray, EdgeExp, coeffs, outarray);
1064  }
1065 
1066 
1067  /** \brief This function evaluates the expansion at a single
1068  * (arbitrary) point of the domain
1069  *
1070  * This function is a wrapper around the virtual function
1071  * \a v_PhysEvaluate()
1072  *
1073  * Based on the value of the expansion at the quadrature
1074  * points provided in \a physvals, this function
1075  * calculates the value of the expansion at an arbitrary
1076  * single points (with coordinates \f$ \mathbf{x_c}\f$
1077  * given by the pointer \a coords). This operation,
1078  * equivalent to \f[ u(\mathbf{x_c}) = \sum_p
1079  * \phi_p(\mathbf{x_c}) \hat{u}_p \f] is evaluated using
1080  * Lagrangian interpolants through the quadrature points:
1081  * \f[ u(\mathbf{x_c}) = \sum_p h_p(\mathbf{x_c}) u_p\f]
1082  *
1083  * \param coords the coordinates of the single point
1084  * \param physvals the interpolated field at the quadrature points
1085  *
1086  * \return returns the value of the expansion at the
1087  * single point
1088  */
1089  NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords,
1090  const Array<OneD, const NekDouble>& physvals)
1091  {
1092  return v_PhysEvaluate(coords,physvals);
1093  }
1094 
1095 
1096  /** \brief This function evaluates the expansion at a single
1097  * (arbitrary) point of the domain
1098  *
1099  * This function is a wrapper around the virtual function
1100  * \a v_PhysEvaluate()
1101  *
1102  * Based on the value of the expansion at the quadrature
1103  * points provided in \a physvals, this function
1104  * calculates the value of the expansion at an arbitrary
1105  * single points associated with the interpolation
1106  * matrices provided in \f$ I \f$.
1107  *
1108  * \param I an Array of lagrange interpolantes evaluated
1109  * at the coordinate and going through the local physical
1110  * quadrature
1111  * \param physvals the interpolated field at the quadrature points
1112  *
1113  * \return returns the value of the expansion at the
1114  * single point
1115  */
1116  NekDouble PhysEvaluate(const Array<OneD, DNekMatSharedPtr>& I,
1117  const Array<OneD, const NekDouble >& physvals)
1118  {
1119  return v_PhysEvaluate(I,physvals);
1120  }
1121 
1122 
1123  /**
1124  * \brief Convert local cartesian coordinate \a xi into local
1125  * collapsed coordinates \a eta
1126  **/
1127  void LocCoordToLocCollapsed(const Array<OneD, const NekDouble>& xi,
1128  Array<OneD, NekDouble>& eta)
1129  {
1130  v_LocCoordToLocCollapsed(xi,eta);
1131  }
1132 
1133  const boost::shared_ptr<SpatialDomains::GeomFactors>& GetMetricInfo(void) const
1134  {
1135  return v_GetMetricInfo();
1136  }
1137 
1138  /// \brief Get the element id of this expansion when used
1139  /// in a list by returning value of #m_elmt_id
1140  STD_REGIONS_EXPORT virtual int v_GetElmtId();
1141 
1142  STD_REGIONS_EXPORT virtual const Array<OneD, const NekDouble>& v_GetPhysNormals(void);
1143 
1144  STD_REGIONS_EXPORT virtual void v_SetPhysNormals(Array<OneD, const NekDouble> &normal);
1145 
1146  STD_REGIONS_EXPORT virtual void v_SetUpPhysNormals(const int edge);
1147 
1148  STD_REGIONS_EXPORT virtual int v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset);
1149 
1150  /**
1151  * @brief Unpack data from input file assuming it comes from the
1152  * same expansion type.
1153  * @see StdExpansion::ExtractDataToCoeffs
1154  */
1155  STD_REGIONS_EXPORT virtual void v_ExtractDataToCoeffs(const NekDouble *data,
1156  const std::vector<unsigned int > &nummodes,
1157  const int nmode_offset,
1158  NekDouble *coeffs);
1159 
1160  STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, Array< OneD, NekDouble> &outarray);
1161 
1162  STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, const Array<OneD, const NekDouble> &Fz, Array< OneD, NekDouble> &outarray);
1163 
1165 
1167 
1168 
1170 
1172 
1174 
1176 
1177  /** \brief Function to evaluate the discrete \f$ L_\infty\f$
1178  * error \f$ |\epsilon|_\infty = \max |u - u_{exact}|\f$ where \f$
1179  * u_{exact}\f$ is given by the array \a sol.
1180  *
1181  * This function takes the physical value space array \a m_phys as
1182  * approximate solution
1183  *
1184  * \param sol array of solution function at physical quadrature
1185  * points
1186  * \return returns the \f$ L_\infty \f$ error as a NekDouble.
1187  */
1188  STD_REGIONS_EXPORT NekDouble Linf(const Array<OneD, const NekDouble>& phys, const Array<OneD, const NekDouble>& sol = NullNekDouble1DArray);
1189 
1190  /** \brief Function to evaluate the discrete \f$ L_2\f$ error,
1191  * \f$ | \epsilon |_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2
1192  * dx \right]^{1/2} d\xi_1 \f$ where \f$ u_{exact}\f$ is given by
1193  * the array \a sol.
1194  *
1195  * This function takes the physical value space array \a m_phys as
1196  * approximate solution
1197  *
1198  * \param sol array of solution function at physical quadrature
1199  * points
1200  * \return returns the \f$ L_2 \f$ error as a double.
1201  */
1202  STD_REGIONS_EXPORT NekDouble L2(const Array<OneD, const NekDouble>& phys, const Array<OneD, const NekDouble>& sol = NullNekDouble1DArray);
1203 
1204  /** \brief Function to evaluate the discrete \f$ H^1\f$
1205  * error, \f$ | \epsilon |^1_{2} = \left [ \int^1_{-1} [u -
1206  * u_{exact}]^2 + \nabla(u - u_{exact})\cdot\nabla(u -
1207  * u_{exact})\cdot dx \right]^{1/2} d\xi_1 \f$ where \f$
1208  * u_{exact}\f$ is given by the array \a sol.
1209  *
1210  * This function takes the physical value space array
1211  * \a m_phys as approximate solution
1212  *
1213  * \param sol array of solution function at physical quadrature
1214  * points
1215  * \return returns the \f$ H_1 \f$ error as a double.
1216  */
1217  STD_REGIONS_EXPORT NekDouble H1(const Array<OneD, const NekDouble>& phys, const Array<OneD, const NekDouble>& sol = NullNekDouble1DArray);
1218 
1219  // I/O routines
1220  const NormalVector & GetEdgeNormal(const int edge) const
1221  {
1222  return v_GetEdgeNormal(edge);
1223  }
1224 
1225  void ComputeEdgeNormal(const int edge)
1226  {
1227  v_ComputeEdgeNormal(edge);
1228  }
1229 
1230  void NegateEdgeNormal(const int edge)
1231  {
1232  v_NegateEdgeNormal(edge);
1233  }
1234 
1235  bool EdgeNormalNegated(const int edge)
1236  {
1237  return v_EdgeNormalNegated(edge);
1238  }
1239 
1240  void ComputeFaceNormal(const int face)
1241  {
1242  v_ComputeFaceNormal(face);
1243  }
1244 
1245  void NegateFaceNormal(const int face)
1246  {
1247  v_NegateFaceNormal(face);
1248  }
1249 
1250  void ComputeVertexNormal(const int vertex)
1251  {
1252  v_ComputeVertexNormal(vertex);
1253  }
1254 
1255  const NormalVector & GetFaceNormal(const int face) const
1256  {
1257  return v_GetFaceNormal(face);
1258  }
1259 
1260  const NormalVector & GetVertexNormal(const int vertex) const
1261  {
1262  return v_GetVertexNormal(vertex);
1263  }
1264 
1265  const NormalVector & GetSurfaceNormal(const int id) const
1266  {
1267  return v_GetSurfaceNormal(id);
1268  }
1269 
1271  {
1273  for (int i = 0; i < m_base.num_elements(); ++i)
1274  {
1275  p.push_back(m_base[i]->GetPointsKey());
1276  }
1277  return p;
1278  }
1279 
1280  STD_REGIONS_EXPORT Array<OneD, unsigned int>
1282  {
1283  return v_GetEdgeInverseBoundaryMap(eid);
1284  }
1285 
1286  STD_REGIONS_EXPORT Array<OneD, unsigned int>
1288  {
1289  return v_GetFaceInverseBoundaryMap(fid,faceOrient);
1290  }
1291 
1293  const DNekScalMatSharedPtr & m_transformationmatrix)
1294  {
1296  m_transformationmatrix);
1297  }
1298 
1299  template<class T>
1300  boost::shared_ptr<T> as()
1301  {
1302 #if defined __INTEL_COMPILER && BOOST_VERSION > 105200
1303  typedef typename boost::shared_ptr<T>::element_type E;
1304  E * p = dynamic_cast< E* >( shared_from_this().get() );
1305  ASSERTL1(p, "Cannot perform cast");
1306  return boost::shared_ptr<T>( shared_from_this(), p );
1307 #else
1308  return boost::dynamic_pointer_cast<T>( shared_from_this() );
1309 #endif
1310  }
1311 
1312  protected:
1313  Array<OneD, LibUtilities::BasisSharedPtr> m_base; /**< Bases needed for the expansion */
1315  int m_ncoeffs; /**< Total number of coefficients used in the expansion */
1319 
1321  {
1322  return v_CreateStdMatrix(mkey);
1323  }
1324 
1325  /** \brief Create the static condensation of a matrix when
1326  using a boundary interior decomposition
1327 
1328  If a matrix system can be represented by
1329  \f$ Mat = \left [ \begin{array}{cc}
1330  A & B \\
1331  C & D \end{array} \right ] \f$
1332  This routine creates a matrix containing the statically
1333  condense system of the form
1334  \f$ Mat = \left [ \begin{array}{cc}
1335  A - B D^{-1} C & B D^{-1} \\
1336  D^{-1} C & D^{-1} \end{array} \right ] \f$
1337  **/
1339 
1340  /** \brief Create an IndexMap which contains mapping information linking any specific
1341  element shape with either its boundaries, edges, faces, verteces, etc.
1342 
1343  The index member of the IndexMapValue struct gives back an integer associated with an entity index
1344  The sign member of the same struct gives back a sign to algebrically apply entities orientation
1345  **/
1347 
1348  STD_REGIONS_EXPORT void BwdTrans_MatOp(const Array<OneD, const NekDouble>& inarray,
1349  Array<OneD, NekDouble> &outarray);
1350 
1351  void BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
1352  Array<OneD, NekDouble> &outarray)
1353  {
1354  v_BwdTrans_SumFac(inarray,outarray);
1355  }
1356 
1357  void IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
1358  Array<OneD, NekDouble> &outarray)
1359  {
1360  v_IProductWRTBase_SumFac(inarray,outarray);
1361  }
1362 
1363  void IProductWRTDerivBase_SumFac(const int dir,
1364  const Array<OneD, const NekDouble>& inarray,
1365  Array<OneD, NekDouble> &outarray)
1366  {
1367  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
1368  }
1369 
1370  // The term _MatFree denotes that the action of the MatrixOperation
1371  // is done withouth actually using the matrix (which then needs to be stored/calculated).
1372  // Although this does not strictly mean that no matrix operations are involved in the
1373  // evaluation of the operation, we use this term in the same context used as in the following
1374  // paper:
1375  // R. C. Kirby, M. G. Knepley, A. Logg, and L. R. Scott,
1376  // "Optimizing the evaluation of finite element matrices," SISC 27:741-758 (2005)
1377  STD_REGIONS_EXPORT void GeneralMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1378  Array<OneD,NekDouble> &outarray,
1379  const StdMatrixKey &mkey);
1380 
1381  STD_REGIONS_EXPORT void MassMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1382  Array<OneD,NekDouble> &outarray,
1383  const StdMatrixKey &mkey);
1384 
1385  void LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1386  Array<OneD,NekDouble> &outarray,
1387  const StdMatrixKey &mkey)
1388  {
1389  v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1390  }
1391 
1393  const Array<OneD, const NekDouble> &inarray,
1394  Array<OneD, NekDouble> &outarray,
1395  Array<OneD, NekDouble> &wsp)
1396  {
1397  v_LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp);
1398  }
1399 
1400  STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree_GenericImpl(const Array<OneD, const NekDouble> &inarray,
1401  Array<OneD,NekDouble> &outarray,
1402  const StdMatrixKey &mkey);
1403 
1404  STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree(const int k1, const int k2,
1405  const Array<OneD, const NekDouble> &inarray,
1406  Array<OneD,NekDouble> &outarray,
1407  const StdMatrixKey &mkey);
1408 
1410  const Array<OneD, const NekDouble> &inarray,
1411  Array<OneD,NekDouble> &outarray,
1412  const StdMatrixKey &mkey);
1413 
1414  STD_REGIONS_EXPORT void WeakDirectionalDerivMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1415  Array<OneD,NekDouble> &outarray,
1416  const StdMatrixKey &mkey);
1417 
1418  STD_REGIONS_EXPORT void MassLevelCurvatureMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1419  Array<OneD,NekDouble> &outarray,
1420  const StdMatrixKey &mkey);
1421 
1422  STD_REGIONS_EXPORT void LinearAdvectionDiffusionReactionMatrixOp_MatFree( const Array<OneD, const NekDouble> &inarray,
1423  Array<OneD,NekDouble> &outarray,
1424  const StdMatrixKey &mkey,
1425  bool addDiffusionTerm = true);
1426 
1427  void HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1428  Array<OneD,NekDouble> &outarray,
1429  const StdMatrixKey &mkey)
1430  {
1431  v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1432  }
1433 
1434  STD_REGIONS_EXPORT void HelmholtzMatrixOp_MatFree_GenericImpl(const Array<OneD, const NekDouble> &inarray,
1435  Array<OneD,NekDouble> &outarray,
1436  const StdMatrixKey &mkey);
1437 
1439  Array<OneD, const NekDouble> &inarray,
1440  Array<OneD, NekDouble> &outarray);
1441 
1443  Array<OneD, NekDouble> &coeffs,
1445 
1447  const Array<OneD, const NekDouble> &Lcoord,
1448  const Array<OneD, const NekDouble> &physvals);
1449 
1450  STD_REGIONS_EXPORT virtual void v_AddEdgeNormBoundaryInt(const int edge,
1451  boost::shared_ptr<StdExpansion> &EdgeExp,
1452  const Array<OneD, const NekDouble> &Fx,
1453  const Array<OneD, const NekDouble> &Fy,
1454  Array<OneD, NekDouble> &outarray);
1455 
1456  STD_REGIONS_EXPORT virtual void v_AddEdgeNormBoundaryInt(const int edge,
1457  boost::shared_ptr<StdExpansion> &EdgeExp,
1458  const Array<OneD, const NekDouble> &Fn,
1459  Array<OneD, NekDouble> &outarray);
1460 
1461  STD_REGIONS_EXPORT virtual void v_AddFaceNormBoundaryInt(const int face,
1462  boost::shared_ptr<StdExpansion> &FaceExp,
1463  const Array<OneD, const NekDouble> &Fn,
1464  Array<OneD, NekDouble> &outarray);
1465 
1466  private:
1467  // Virtual functions
1468  STD_REGIONS_EXPORT virtual int v_GetNverts() const = 0;
1469  STD_REGIONS_EXPORT virtual int v_GetNedges() const;
1470 
1471  STD_REGIONS_EXPORT virtual int v_GetNfaces() const;
1472 
1473  STD_REGIONS_EXPORT virtual int v_NumBndryCoeffs() const;
1474 
1475  STD_REGIONS_EXPORT virtual int v_NumDGBndryCoeffs() const;
1476 
1477  STD_REGIONS_EXPORT virtual int v_GetEdgeNcoeffs(const int i) const;
1478 
1479  STD_REGIONS_EXPORT virtual int v_GetTotalEdgeIntNcoeffs() const;
1480 
1481  STD_REGIONS_EXPORT virtual int v_GetEdgeNumPoints(const int i) const;
1482 
1483  STD_REGIONS_EXPORT virtual int v_DetCartesianDirOfEdge(const int edge);
1484 
1485  STD_REGIONS_EXPORT virtual const LibUtilities::BasisKey v_DetEdgeBasisKey(const int i) const;
1486 
1487  STD_REGIONS_EXPORT virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const;
1488 
1489  STD_REGIONS_EXPORT virtual int v_GetFaceNumPoints(const int i) const;
1490 
1491  STD_REGIONS_EXPORT virtual int v_GetFaceNcoeffs(const int i) const;
1492 
1493  STD_REGIONS_EXPORT virtual int v_GetFaceIntNcoeffs(const int i) const;
1494 
1495  STD_REGIONS_EXPORT virtual int v_GetTotalFaceIntNcoeffs() const;
1496 
1497  STD_REGIONS_EXPORT virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const;
1498 
1500 
1502 
1503  STD_REGIONS_EXPORT virtual int v_GetShapeDimension() const;
1504 
1506 
1507  STD_REGIONS_EXPORT virtual void v_BwdTrans (const Array<OneD, const NekDouble>& inarray,
1508  Array<OneD, NekDouble> &outarray) = 0;
1509 
1510  /**
1511  * @brief Transform a given function from physical quadrature space
1512  * to coefficient space.
1513  * @see StdExpansion::FwdTrans
1514  */
1515  STD_REGIONS_EXPORT virtual void v_FwdTrans (
1516  const Array<OneD, const NekDouble>& inarray,
1517  Array<OneD, NekDouble> &outarray) = 0;
1518 
1519  /**
1520  * @brief Calculates the inner product of a given function \a f
1521  * with the different modes of the expansion
1522  */
1524  const Array<OneD, const NekDouble>& inarray,
1525  Array<OneD, NekDouble> &outarray) = 0;
1526 
1528  const Array<OneD, const NekDouble>& base,
1529  const Array<OneD, const NekDouble>& inarray,
1530  Array<OneD, NekDouble>& outarray,
1531  int coll_check)
1532  {
1533  ASSERTL0(false, "StdExpansion::v_IProductWRTBase has no (and should have no) implementation");
1534  }
1535 
1536  STD_REGIONS_EXPORT virtual void v_IProductWRTDerivBase (const int dir,
1537  const Array<OneD, const NekDouble>& inarray,
1538  Array<OneD, NekDouble> &outarray);
1539 
1540  STD_REGIONS_EXPORT virtual void v_FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
1541  Array<OneD, NekDouble> &outarray);
1542 
1543  STD_REGIONS_EXPORT virtual NekDouble v_Integral(const Array<OneD, const NekDouble>& inarray );
1544 
1545  STD_REGIONS_EXPORT virtual void v_PhysDeriv (const Array<OneD, const NekDouble>& inarray,
1546  Array<OneD, NekDouble> &out_d1,
1547  Array<OneD, NekDouble> &out_d2,
1548  Array<OneD, NekDouble> &out_d3);
1549 
1550  STD_REGIONS_EXPORT virtual void v_PhysDeriv_s (const Array<OneD, const NekDouble>& inarray,
1551  Array<OneD, NekDouble> &out_ds);
1552 
1553  STD_REGIONS_EXPORT virtual void v_PhysDeriv_n(const Array<OneD, const NekDouble>& inarray,
1554  Array<OneD, NekDouble>& out_dn);
1555  STD_REGIONS_EXPORT virtual void v_PhysDeriv(const int dir,
1556  const Array<OneD, const NekDouble>& inarray,
1557  Array<OneD, NekDouble> &out_d0);
1558 
1559  STD_REGIONS_EXPORT virtual void v_PhysDirectionalDeriv(const Array<OneD, const NekDouble>& inarray,
1560  const Array<OneD, const NekDouble>& direction,
1561  Array<OneD, NekDouble> &outarray);
1562 
1563  STD_REGIONS_EXPORT virtual void v_StdPhysDeriv (const Array<OneD, const NekDouble>& inarray,
1564  Array<OneD, NekDouble> &out_d1,
1565  Array<OneD, NekDouble> &out_d2,
1566  Array<OneD, NekDouble> &out_d3);
1567 
1568  STD_REGIONS_EXPORT virtual void v_StdPhysDeriv (const int dir,
1569  const Array<OneD, const NekDouble>& inarray,
1570  Array<OneD, NekDouble> &outarray);
1571 
1572  STD_REGIONS_EXPORT virtual void v_AddRobinMassMatrix(const int edgeid, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat);
1573 
1574  STD_REGIONS_EXPORT virtual void v_AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble> &primCoeffs, Array<OneD, NekDouble> &coeffs);
1575 
1576  STD_REGIONS_EXPORT virtual void v_DGDeriv(
1577  const int dir,
1578  const Array<OneD, const NekDouble>& inarray,
1579  Array<OneD, boost::shared_ptr<StdExpansion> > &EdgeExp,
1580  Array<OneD, Array<OneD, NekDouble> > &coeffs,
1581  Array<OneD, NekDouble> &outarray);
1582 
1583  STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
1584 
1585  STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, DNekMatSharedPtr >& I, const Array<OneD, const NekDouble> & physvals);
1586 
1588  const Array<OneD, const NekDouble>& xi,
1589  Array<OneD, NekDouble>& eta);
1590 
1591 
1592  STD_REGIONS_EXPORT virtual void v_FillMode(const int mode, Array<OneD, NekDouble> &outarray);
1593 
1595 
1597 
1598  STD_REGIONS_EXPORT virtual void v_GetCoords(Array<OneD, NekDouble> &coords_0,
1599  Array<OneD, NekDouble> &coords_1,
1600  Array<OneD, NekDouble> &coords_2);
1601 
1602  STD_REGIONS_EXPORT virtual void v_GetCoord(const Array<OneD, const NekDouble>& Lcoord,
1603  Array<OneD, NekDouble> &coord);
1604 
1605  STD_REGIONS_EXPORT virtual int v_GetCoordim(void);
1606 
1607  STD_REGIONS_EXPORT virtual void v_GetBoundaryMap(Array<OneD, unsigned int>& outarray);
1608 
1609  STD_REGIONS_EXPORT virtual void v_GetInteriorMap(Array<OneD, unsigned int>& outarray);
1610 
1611  STD_REGIONS_EXPORT virtual int v_GetVertexMap(int localVertexId,
1612  bool useCoeffPacking = false);
1613 
1614  STD_REGIONS_EXPORT virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient,
1615  Array<OneD, unsigned int> &maparray,
1616  Array<OneD, int> &signarray);
1617 
1618  STD_REGIONS_EXPORT virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient,
1619  Array<OneD, unsigned int> &maparray,
1620  Array<OneD, int> &signarray);
1621 
1622  STD_REGIONS_EXPORT virtual void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient,
1623  Array<OneD, unsigned int> &maparray,
1624  Array<OneD, int> &signarray);
1625 
1626  STD_REGIONS_EXPORT virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient,
1627  Array<OneD, unsigned int> &maparray,
1628  Array<OneD, int> &signarray,
1629  int nummodesA = -1, int nummodesB = -1);
1630 
1631  /**
1632  * @brief Extract the physical values along edge \a edge from \a
1633  * inarray into \a outarray following the local edge orientation
1634  * and point distribution defined by defined in \a EdgeExp.
1635  */
1636  STD_REGIONS_EXPORT virtual void v_GetEdgePhysVals(const int edge, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray);
1637 
1638  STD_REGIONS_EXPORT virtual void v_GetEdgePhysVals(const int edge, const boost::shared_ptr<StdExpansion> &EdgeExp, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray);
1639 
1640  STD_REGIONS_EXPORT virtual void v_GetTracePhysVals(const int edge, const boost::shared_ptr<StdExpansion> &EdgeExp, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray, StdRegions::Orientation orient = eNoOrientation);
1641 
1642  STD_REGIONS_EXPORT virtual void v_GetVertexPhysVals(const int vertex, const Array<OneD, const NekDouble> &inarray, NekDouble &outarray);
1643 
1644  STD_REGIONS_EXPORT virtual void v_GetEdgeInterpVals(const int edge,
1645  const Array<OneD, const NekDouble> &inarray,Array<OneD,NekDouble> &outarray);
1646 
1648  const int edge,
1649  Array<OneD, NekDouble> &outarray);
1650 
1652  const int face,
1653  const boost::shared_ptr<StdExpansion> &FaceExp,
1654  const Array<OneD, const NekDouble> &inarray,
1655  Array<OneD, NekDouble> &outarray,
1656  StdRegions::Orientation orient);
1657 
1659  const Array<OneD, const NekDouble> &inarray,
1660  Array<OneD, NekDouble> &outarray);
1661 
1663  const Array<OneD, const NekDouble> &inarray,
1664  Array<OneD, NekDouble> &outarray);
1665 
1666  STD_REGIONS_EXPORT virtual const boost::shared_ptr<SpatialDomains::GeomFactors>& v_GetMetricInfo() const;
1667 
1668  STD_REGIONS_EXPORT virtual void v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
1669  Array<OneD, NekDouble> &outarray);
1670 
1671  STD_REGIONS_EXPORT virtual void v_IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
1672  Array<OneD, NekDouble> &outarray);
1673 
1674  STD_REGIONS_EXPORT virtual void v_IProductWRTDerivBase_SumFac(const int dir,
1675  const Array<OneD, const NekDouble>& inarray,
1676  Array<OneD, NekDouble> &outarray);
1677 
1678  STD_REGIONS_EXPORT virtual void v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
1679  Array<OneD,NekDouble> &outarray,
1680  const StdMatrixKey &mkey);
1681 
1682  STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
1683  Array<OneD,NekDouble> &outarray,
1684  const StdMatrixKey &mkey);
1685 
1686  STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD,NekDouble> &array,
1687  const StdMatrixKey &mkey);
1688 
1690  int numMin,
1691  const Array<OneD, const NekDouble> &inarray,
1692  Array<OneD,NekDouble> &outarray);
1693 
1694  STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(const int k1, const int k2,
1695  const Array<OneD, const NekDouble> &inarray,
1696  Array<OneD,NekDouble> &outarray,
1697  const StdMatrixKey &mkey);
1698 
1699  STD_REGIONS_EXPORT virtual void v_WeakDerivMatrixOp(const int i,
1700  const Array<OneD, const NekDouble> &inarray,
1701  Array<OneD,NekDouble> &outarray,
1702  const StdMatrixKey &mkey);
1703 
1704  STD_REGIONS_EXPORT virtual void v_WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
1705  Array<OneD,NekDouble> &outarray,
1706  const StdMatrixKey &mkey);
1707 
1708  STD_REGIONS_EXPORT virtual void v_MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
1709  Array<OneD,NekDouble> &outarray,
1710  const StdMatrixKey &mkey);
1711 
1713  const NekDouble> &inarray,
1714  Array<OneD,NekDouble> &outarray,
1715  const StdMatrixKey &mkey,
1716  bool addDiffusionTerm=true);
1717 
1718  STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
1719  Array<OneD,NekDouble> &outarray,
1720  const StdMatrixKey &mkey);
1721 
1722  STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1723  Array<OneD,NekDouble> &outarray,
1724  const StdMatrixKey &mkey);
1725 
1727  const Array<OneD, const NekDouble> &inarray,
1728  Array<OneD, NekDouble> &outarray,
1729  Array<OneD, NekDouble> &wsp);
1730 
1731  STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1732  Array<OneD,NekDouble> &outarray,
1733  const StdMatrixKey &mkey);
1734 
1735  STD_REGIONS_EXPORT virtual const NormalVector & v_GetEdgeNormal(const int edge) const;
1736 
1737  STD_REGIONS_EXPORT virtual void v_ComputeEdgeNormal(const int edge);
1738 
1739  STD_REGIONS_EXPORT virtual void v_NegateEdgeNormal(const int edge);
1740 
1741  STD_REGIONS_EXPORT virtual bool v_EdgeNormalNegated(const int edge);
1742 
1743  STD_REGIONS_EXPORT virtual void v_ComputeFaceNormal(const int face);
1744 
1745  STD_REGIONS_EXPORT virtual void v_NegateFaceNormal(const int face);
1746 
1747  STD_REGIONS_EXPORT virtual const NormalVector & v_GetVertexNormal(const int vertex) const;
1748 
1749  STD_REGIONS_EXPORT virtual void v_ComputeVertexNormal(const int vertex);
1750 
1751  STD_REGIONS_EXPORT virtual const NormalVector & v_GetFaceNormal(const int face) const;
1752  STD_REGIONS_EXPORT virtual const NormalVector &
1753  v_GetSurfaceNormal(const int id) const;
1754 
1755  STD_REGIONS_EXPORT virtual Array<OneD, unsigned int>
1756  v_GetEdgeInverseBoundaryMap(int eid);
1757 
1758  STD_REGIONS_EXPORT virtual Array<OneD, unsigned int>
1760 
1762  const DNekScalMatSharedPtr & m_transformationmatrix);
1763 
1764  };
1765 
1766 
1767  typedef boost::shared_ptr<StdExpansion> StdExpansionSharedPtr;
1768  typedef std::vector< StdExpansionSharedPtr > StdExpansionVector;
1770 
1771  /**
1772  * This function is a wrapper around the virtual function
1773  * \a v_FwdTrans()
1774  *
1775  * Given a function evaluated at the quadrature points, this
1776  * function calculates the expansion coefficients such that the
1777  * resulting expansion approximates the original function.
1778  *
1779  * The calculation of the expansion coefficients is done using a
1780  * Galerkin projection. This is equivalent to the operation:
1781  * \f[ \mathbf{\hat{u}} = \mathbf{M}^{-1} \mathbf{I}\f]
1782  * where
1783  * - \f$\mathbf{M}[p][q]= \int\phi_p(\mathbf{\xi})\phi_q(
1784  * \mathbf{\xi}) d\mathbf{\xi}\f$ is the Mass matrix
1785  * - \f$\mathbf{I}[p] = \int\phi_p(\mathbf{\xi}) u(\mathbf{\xi})
1786  * d\mathbf{\xi}\f$
1787  *
1788  * This function takes the array \a inarray as the values of the
1789  * function evaluated at the quadrature points
1790  * (i.e. \f$\mathbf{u}\f$),
1791  * and stores the resulting coefficients \f$\mathbf{\hat{u}}\f$
1792  * in the \a outarray
1793  *
1794  * @param inarray array of the function discretely evaluated at the
1795  * quadrature points
1796  *
1797  * @param outarray array of the function coefficieints
1798  */
1799  inline void StdExpansion::FwdTrans (const Array<OneD, const NekDouble>& inarray,
1800  Array<OneD, NekDouble> &outarray)
1801  {
1802  v_FwdTrans(inarray,outarray);
1803  }
1804 
1805  } //end of namespace
1806 } //end of namespace
1807 
1808 #endif //STANDARDDEXPANSION_H