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  /**
417  * @brief Returns the number of trace elements connected to this
418  * element.
419  *
420  * For example, a quadrilateral has four edges, so this function
421  * would return 4.
422  */
423  int GetNtrace() const
424  {
425  const int nBase = m_base.num_elements();
426  return
427  nBase == 1 ? 2 :
428  nBase == 2 ? GetNedges() :
429  nBase == 3 ? GetNfaces() : 0;
430  }
431 
432  /** \brief This function returns the shape of the expansion domain
433  *
434  * This function is a wrapper around the virtual function
435  * \a v_DetShapeType()
436  *
437  * The different shape types implemented in the code are defined
438  * in the ::ShapeType enumeration list. As a result, the
439  * function will return one of the types of this enumeration list.
440  *
441  * \return returns the shape of the expansion domain
442  */
444  {
445  return v_DetShapeType();
446  }
447 
448  int GetShapeDimension() const
449  {
450  return v_GetShapeDimension();
451  }
452 
454  {
456  }
457 
458 
459  /** \brief This function performs the Backward transformation from
460  * coefficient space to physical space
461  *
462  * This function is a wrapper around the virtual function
463  * \a v_BwdTrans()
464  *
465  * Based on the expansion coefficients, this function evaluates the
466  * expansion at the quadrature points. This is equivalent to the
467  * operation \f[ u(\xi_{1i}) =
468  * \sum_{p=0}^{P-1} \hat{u}_p \phi_p(\xi_{1i}) \f] which can be
469  * evaluated as \f$ {\bf u} = {\bf B}^T {\bf \hat{u}} \f$ with
470  * \f${\bf B}[i][j] = \phi_i(\xi_{j})\f$
471  *
472  * This function requires that the coefficient array
473  * \f$\mathbf{\hat{u}}\f$ provided as \a inarray.
474  *
475  * The resulting array
476  * \f$\mathbf{u}[m]=u(\mathbf{\xi}_m)\f$ containing the
477  * expansion evaluated at the quadrature points, is stored
478  * in the \a outarray.
479  *
480  * \param inarray contains the values of the expansion
481  * coefficients (input of the function)
482  *
483  * \param outarray contains the values of the expansion evaluated
484  * at the quadrature points (output of the function)
485  */
486 
487  void BwdTrans (const Array<OneD, const NekDouble>& inarray,
488  Array<OneD, NekDouble> &outarray)
489  {
490  v_BwdTrans (inarray, outarray);
491  }
492 
493  /**
494  * @brief This function performs the Forward transformation from
495  * physical space to coefficient space.
496  */
497  inline void FwdTrans (const Array<OneD, const NekDouble>& inarray,
498  Array<OneD, NekDouble> &outarray);
499 
500  void FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
501  Array<OneD, NekDouble> &outarray)
502  {
503  v_FwdTrans_BndConstrained(inarray,outarray);
504  }
505 
506  /** \brief This function integrates the specified function over the
507  * domain
508  *
509  * This function is a wrapper around the virtual function
510  * \a v_Integral()
511  *
512  * Based on the values of the function evaluated at the quadrature
513  * points (which are stored in \a inarray), this function calculates
514  * the integral of this function over the domain. This is
515  * equivalent to the numerical evaluation of the operation
516  * \f[ I=\int u(\mathbf{\xi})d \mathbf{\xi}\f]
517  *
518  * \param inarray values of the function to be integrated evaluated
519  * at the quadrature points (i.e.
520  * \a inarray[m]=\f$u(\mathbf{\xi}_m)\f$)
521  * \return returns the value of the calculated integral
522  *
523  * Inputs:\n
524 
525  - \a inarray: definition of function to be returned at quadrature point
526  of expansion.
527 
528  Outputs:\n
529 
530  - returns \f$\int^1_{-1}\int^1_{-1} u(\xi_1, \xi_2) J[i,j] d
531  \xi_1 d \xi_2 \f$ where \f$inarray[i,j] =
532  u(\xi_{1i},\xi_{2j}) \f$ and \f$ J[i,j] \f$ is the
533  Jacobian evaluated at the quadrature point.
534  *
535  */
536  NekDouble Integral(const Array<OneD, const NekDouble>& inarray )
537  {
538  return v_Integral(inarray);
539  }
540 
541  /** \brief This function fills the array \a outarray with the
542  * \a mode-th mode of the expansion
543  *
544  * This function is a wrapper around the virtual function
545  * \a v_FillMode()
546  *
547  * The requested mode is evaluated at the quadrature points
548  *
549  * \param mode the mode that should be filled
550  * \param outarray contains the values of the \a mode-th mode of the
551  * expansion evaluated at the quadrature points (output of the
552  * function)
553  */
554  void FillMode(const int mode, Array<OneD, NekDouble> &outarray)
555  {
556  v_FillMode(mode, outarray);
557  }
558 
559  /** \brief this function calculates the inner product of a given
560  * function \a f with the different modes of the expansion
561  *
562  * This function is a wrapper around the virtual function
563  * \a v_IProductWRTBase()
564  *
565  * This is equivalent to the numerical evaluation of
566  * \f[ I[p] = \int \phi_p(\mathbf{x}) f(\mathbf{x}) d\mathbf{x}\f]
567  * \f$ \begin{array}{rcl} I_{pq} = (\phi_q \phi_q, u) & = &
568  \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \phi_p(\xi_{0,i})
569  \phi_q(\xi_{1,j}) w^0_i w^1_j u(\xi_{0,i} \xi_{1,j})
570  J_{i,j}\\ & = & \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i})
571  \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j}
572  J_{i,j} \end{array} \f$
573 
574  where
575 
576  \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
577 
578  which can be implemented as
579 
580  \f$ f_{qi} = \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j} =
581  {\bf B_1 U} \f$
582  \f$ I_{pq} = \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i}) f_{qi} =
583  {\bf B_0 F} \f$
584  *
585  * \param inarray contains the values of the function \a f
586  * evaluated at the quadrature points
587  * \param outarray contains the values of the inner product of \a f
588  * with the different modes, i.e. \f$ outarray[p] = I[p]\f$
589  * (output of the function)
590  */
591  void IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
592  Array<OneD, NekDouble> &outarray)
593  {
594  v_IProductWRTBase(inarray, outarray);
595  }
596 
598  const Array<OneD, const NekDouble>& base,
599  const Array<OneD, const NekDouble>& inarray,
600  Array<OneD, NekDouble> &outarray,
601  int coll_check)
602  {
603  v_IProductWRTBase(base, inarray, outarray, coll_check);
604  }
605 
606 
607  void IProductWRTDerivBase(const int dir,
608  const Array<OneD, const NekDouble>& inarray,
609  Array<OneD, NekDouble> &outarray)
610  {
611  v_IProductWRTDerivBase(dir,inarray, outarray);
612  }
613 
614  /// \brief Get the element id of this expansion when used
615  /// in a list by returning value of #m_elmt_id
616  inline int GetElmtId()
617  {
618  return m_elmt_id;
619  }
620 
621 
622  /// \brief Set the element id of this expansion when used
623  /// in a list by returning value of #m_elmt_id
624  inline void SetElmtId(const int id)
625  {
626  m_elmt_id = id;
627  }
628 
629  /** \brief this function returns the physical coordinates of the
630  * quadrature points of the expansion
631  *
632  * This function is a wrapper around the virtual function
633  * \a v_GetCoords()
634  *
635  * \param coords an array containing the coordinates of the
636  * quadrature points (output of the function)
637  */
638  void GetCoords(Array<OneD, NekDouble> &coords_1,
639  Array<OneD, NekDouble> &coords_2 = NullNekDouble1DArray,
640  Array<OneD, NekDouble> &coords_3 = NullNekDouble1DArray)
641  {
642  v_GetCoords(coords_1,coords_2,coords_3);
643  }
644 
645  /** \brief given the coordinates of a point of the element in the
646  * local collapsed coordinate system, this function calculates the
647  * physical coordinates of the point
648  *
649  * This function is a wrapper around the virtual function
650  * \a v_GetCoord()
651  *
652  * \param Lcoords the coordinates in the local collapsed
653  * coordinate system
654  * \param coords the physical coordinates (output of the function)
655  */
656  void GetCoord(const Array<OneD, const NekDouble>& Lcoord,
657  Array<OneD, NekDouble> &coord)
658  {
659  v_GetCoord(Lcoord, coord);
660  }
661 
663  {
664  return m_stdMatrixManager[mkey];
665  }
666 
668  {
669  return m_stdStaticCondMatrixManager[mkey];
670  }
671 
673  {
674  return m_IndexMapManager[ikey];
675  }
676 
677  const Array<OneD, const NekDouble>& GetPhysNormals(void)
678  {
679  return v_GetPhysNormals();
680  }
681 
682  void SetPhysNormals(Array<OneD, const NekDouble> &normal)
683  {
684  v_SetPhysNormals(normal);
685  }
686 
687  STD_REGIONS_EXPORT virtual void SetUpPhysNormals(const int edge);
688 
689  void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, Array< OneD, NekDouble> &outarray)
690  {
691  v_NormVectorIProductWRTBase(Fx,Fy,outarray);
692  }
693 
694  void NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, const Array<OneD, const NekDouble> &Fz, Array< OneD, NekDouble> &outarray)
695  {
696  v_NormVectorIProductWRTBase(Fx,Fy,Fz,outarray);
697  }
698 
700  {
701  return v_GetLocStaticCondMatrix(mkey);
702  }
703 
705  {
706  return v_DropLocStaticCondMatrix(mkey);
707  }
708 
710  {
711  return v_GetForient(face);
712  }
713 
715  {
716  return v_GetEorient(edge);
717  }
718 
720  {
721  return v_GetPorient(point);
722  }
723 
725  {
726  return v_GetCartesianEorient(edge);
727  }
728 
730  Array<OneD, NekDouble> &coeffs,
732  {
733  v_SetCoeffsToOrientation(coeffs, dir);
734  }
735 
738  Array<OneD, const NekDouble> &inarray,
739  Array<OneD, NekDouble> &outarray)
740  {
741  v_SetCoeffsToOrientation(dir,inarray,outarray);
742  }
743 
744  int CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset)
745  {
746  return v_CalcNumberOfCoefficients(nummodes,modes_offset);
747  }
748 
749  void ExtractDataToCoeffs(const NekDouble *data,
750  const std::vector<unsigned int > &nummodes,
751  const int nmodes_offset,
752  NekDouble *coeffs)
753  {
754  v_ExtractDataToCoeffs(data,nummodes,nmodes_offset,coeffs);
755  }
756 
757  // virtual functions related to LocalRegions
759  const Array<OneD, const NekDouble> &Lcoord,
760  const Array<OneD, const NekDouble> &physvals);
761 
762 
764  {
765  return v_GetCoordim();
766  }
767 
768  void GetBoundaryMap(Array<OneD, unsigned int> &outarray)
769  {
770  v_GetBoundaryMap(outarray);
771  }
772 
773  void GetInteriorMap(Array<OneD, unsigned int> &outarray)
774  {
775  v_GetInteriorMap(outarray);
776  }
777 
778  int GetVertexMap(const int localVertexId,
779  bool useCoeffPacking = false)
780  {
781  return v_GetVertexMap(localVertexId,useCoeffPacking);
782  }
783 
784  void GetEdgeInteriorMap(const int eid, const Orientation edgeOrient,
785  Array<OneD, unsigned int> &maparray,
786  Array<OneD, int> &signarray)
787  {
788  v_GetEdgeInteriorMap(eid,edgeOrient,maparray,signarray);
789  }
790 
791  void GetFaceInteriorMap(const int fid, const Orientation faceOrient,
792  Array<OneD, unsigned int> &maparray,
793  Array<OneD, int> &signarray)
794  {
795  v_GetFaceInteriorMap(fid,faceOrient,maparray,signarray);
796  }
797 
798  void GetEdgeToElementMap(const int eid, const Orientation edgeOrient,
799  Array<OneD, unsigned int> &maparray,
800  Array<OneD, int> &signarray)
801  {
802  v_GetEdgeToElementMap(eid,edgeOrient,maparray,signarray);
803  }
804 
805  void GetFaceToElementMap(const int fid, const Orientation faceOrient,
806  Array<OneD, unsigned int> &maparray,
807  Array<OneD, int> &signarray,
808  int nummodesA = -1, int nummodesB = -1)
809  {
810  v_GetFaceToElementMap(fid,faceOrient,maparray,signarray,
811  nummodesA,nummodesB);
812  }
813 
814 
815  /**
816  * @brief Extract the physical values along edge \a edge from \a
817  * inarray into \a outarray following the local edge orientation
818  * and point distribution defined by defined in \a EdgeExp.
819  */
820 
821  void GetEdgePhysVals(const int edge, const Array<OneD,
822  const NekDouble> &inarray,
823  Array<OneD,NekDouble> &outarray)
824  {
825  v_GetEdgePhysVals(edge,inarray,outarray);
826  }
827 
828  void GetEdgePhysVals(const int edge,
829  const boost::shared_ptr<StdExpansion> &EdgeExp,
830  const Array<OneD, const NekDouble> &inarray,
831  Array<OneD,NekDouble> &outarray)
832  {
833  v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
834  }
835 
836  void GetTracePhysVals(const int edge, const boost::shared_ptr<StdExpansion> &EdgeExp, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
837  {
838  v_GetTracePhysVals(edge,EdgeExp,inarray,outarray);
839  }
840 
841  void GetVertexPhysVals(const int vertex,
842  const Array<OneD, const NekDouble> &inarray,
843  NekDouble &outarray)
844  {
845  v_GetVertexPhysVals(vertex, inarray, outarray);
846  }
847 
848  void GetEdgeInterpVals(const int edge,const Array<OneD,
849  const NekDouble> &inarray,
850  Array<OneD,NekDouble> &outarray)
851  {
852  v_GetEdgeInterpVals(edge, inarray, outarray);
853  }
854 
855  /**
856  * @brief Extract the metric factors to compute the contravariant
857  * fluxes along edge \a edge and stores them into \a outarray
858  * following the local edge orientation (i.e. anticlockwise
859  * convention).
860  */
862  const int edge,
863  Array<OneD, NekDouble> &outarray)
864  {
865  v_GetEdgeQFactors(edge, outarray);
866  }
867 
868 
869 
871  const int face,
872  const boost::shared_ptr<StdExpansion> &FaceExp,
873  const Array<OneD, const NekDouble> &inarray,
874  Array<OneD, NekDouble> &outarray,
876  {
877  v_GetFacePhysVals(face, FaceExp, inarray, outarray, orient);
878  }
879 
881  const Array<OneD, const NekDouble> &inarray,
882  Array<OneD, NekDouble> &outarray)
883  {
884  v_MultiplyByQuadratureMetric(inarray, outarray);
885  }
886 
888  const Array<OneD, const NekDouble> &inarray,
889  Array<OneD, NekDouble> & outarray)
890  {
891  v_MultiplyByStdQuadratureMetric(inarray, outarray);
892  }
893 
894  // Matrix Routines
895 
896  /** \brief this function generates the mass matrix
897  * \f$\mathbf{M}[i][j] =
898  * \int \phi_i(\mathbf{x}) \phi_j(\mathbf{x}) d\mathbf{x}\f$
899  *
900  * \return returns the mass matrix
901  */
902 
904 
905  STD_REGIONS_EXPORT void GeneralMatrixOp(const Array<OneD, const NekDouble> &inarray,
906  Array<OneD,NekDouble> &outarray,
907  const StdMatrixKey &mkey);
908 
909  void MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
910  Array<OneD,NekDouble> &outarray,
911  const StdMatrixKey &mkey)
912  {
913  v_MassMatrixOp(inarray,outarray,mkey);
914  }
915 
916  void LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
917  Array<OneD,NekDouble> &outarray,
918  const StdMatrixKey &mkey)
919  {
920  v_LaplacianMatrixOp(inarray,outarray,mkey);
921  }
922 
923  void ReduceOrderCoeffs(int numMin,
924  const Array<OneD, const NekDouble> &inarray,
925  Array<OneD,NekDouble> &outarray)
926  {
927  v_ReduceOrderCoeffs(numMin,inarray,outarray);
928  }
929 
930  void SVVLaplacianFilter(Array<OneD,NekDouble> &array,
931  const StdMatrixKey &mkey)
932  {
933  v_SVVLaplacianFilter(array,mkey);
934  }
935 
936  void LaplacianMatrixOp(const int k1, const int k2,
937  const Array<OneD, const NekDouble> &inarray,
938  Array<OneD,NekDouble> &outarray,
939  const StdMatrixKey &mkey)
940  {
941  v_LaplacianMatrixOp(k1,k2,inarray,outarray,mkey);
942  }
943 
944  void WeakDerivMatrixOp(const int i,
945  const Array<OneD, const NekDouble> &inarray,
946  Array<OneD,NekDouble> &outarray,
947  const StdMatrixKey &mkey)
948  {
949  v_WeakDerivMatrixOp(i,inarray,outarray,mkey);
950  }
951 
952  void WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
953  Array<OneD,NekDouble> &outarray,
954  const StdMatrixKey &mkey)
955  {
956  v_WeakDirectionalDerivMatrixOp(inarray,outarray,mkey);
957  }
958 
959  void MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
960  Array<OneD,NekDouble> &outarray,
961  const StdMatrixKey &mkey)
962  {
963  v_MassLevelCurvatureMatrixOp(inarray,outarray,mkey);
964  }
965 
966  void LinearAdvectionDiffusionReactionMatrixOp(const Array<OneD, const NekDouble> &inarray,
967  Array<OneD,NekDouble> &outarray,
968  const StdMatrixKey &mkey,
969  bool addDiffusionTerm = true)
970  {
971  v_LinearAdvectionDiffusionReactionMatrixOp(inarray,outarray,mkey,addDiffusionTerm);
972  }
973 
974  /**
975  * @param inarray Input array @f$ \mathbf{u} @f$.
976  * @param outarray Output array @f$ \boldsymbol{\nabla^2u}
977  * + \lambda \boldsymbol{u} @f$.
978  * @param mkey
979  */
980  void HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
981  Array<OneD,NekDouble> &outarray,
982  const StdMatrixKey &mkey)
983  {
984  v_HelmholtzMatrixOp(inarray,outarray,mkey);
985  }
986 
988  {
989  return v_GenMatrix(mkey);
990  }
991 
992  void PhysDeriv (const Array<OneD, const NekDouble>& inarray,
993  Array<OneD, NekDouble> &out_d0,
994  Array<OneD, NekDouble> &out_d1 = NullNekDouble1DArray,
995  Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray)
996  {
997  v_PhysDeriv (inarray, out_d0, out_d1, out_d2);
998  }
999 
1000  void PhysDeriv(const int dir,
1001  const Array<OneD, const NekDouble>& inarray,
1002  Array<OneD, NekDouble> &outarray)
1003  {
1004  v_PhysDeriv (dir, inarray, outarray);
1005  }
1006 
1007  void PhysDeriv_s(const Array<OneD, const NekDouble>& inarray,
1008  Array<OneD, NekDouble> &out_ds)
1009  {
1010  v_PhysDeriv_s(inarray,out_ds);
1011  }
1012 
1013  void PhysDeriv_n(const Array<OneD, const NekDouble>& inarray,
1014  Array<OneD, NekDouble>& out_dn)
1015  {
1016  v_PhysDeriv_n(inarray,out_dn);
1017  }
1018 
1019  void PhysDirectionalDeriv(const Array<OneD, const NekDouble>& inarray,
1020  const Array<OneD, const NekDouble>& direction,
1021  Array<OneD, NekDouble> &outarray)
1022  {
1023  v_PhysDirectionalDeriv (inarray, direction, outarray);
1024  }
1025 
1026  void StdPhysDeriv(const Array<OneD, const NekDouble>& inarray,
1027  Array<OneD, NekDouble> &out_d0,
1028  Array<OneD, NekDouble> &out_d1 = NullNekDouble1DArray,
1029  Array<OneD, NekDouble> &out_d2 = NullNekDouble1DArray)
1030  {
1031  v_StdPhysDeriv(inarray, out_d0, out_d1, out_d2);
1032  }
1033 
1034  void StdPhysDeriv (const int dir,
1035  const Array<OneD, const NekDouble>& inarray,
1036  Array<OneD, NekDouble> &outarray)
1037  {
1038  v_StdPhysDeriv(dir,inarray,outarray);
1039  }
1040 
1041  void AddRobinMassMatrix(const int edgeid, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
1042  {
1043  v_AddRobinMassMatrix(edgeid,primCoeffs,inoutmat);
1044  }
1045 
1046  void AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble> &primCoeffs, Array<OneD, NekDouble> &coeffs)
1047  {
1048  v_AddRobinEdgeContribution(edgeid, primCoeffs, coeffs);
1049  }
1050 
1051  /** \brief This function evaluates the expansion at a single
1052  * (arbitrary) point of the domain
1053  *
1054  * This function is a wrapper around the virtual function
1055  * \a v_PhysEvaluate()
1056  *
1057  * Based on the value of the expansion at the quadrature
1058  * points provided in \a physvals, this function
1059  * calculates the value of the expansion at an arbitrary
1060  * single points (with coordinates \f$ \mathbf{x_c}\f$
1061  * given by the pointer \a coords). This operation,
1062  * equivalent to \f[ u(\mathbf{x_c}) = \sum_p
1063  * \phi_p(\mathbf{x_c}) \hat{u}_p \f] is evaluated using
1064  * Lagrangian interpolants through the quadrature points:
1065  * \f[ u(\mathbf{x_c}) = \sum_p h_p(\mathbf{x_c}) u_p\f]
1066  *
1067  * \param coords the coordinates of the single point
1068  * \param physvals the interpolated field at the quadrature points
1069  *
1070  * \return returns the value of the expansion at the
1071  * single point
1072  */
1073  NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords,
1074  const Array<OneD, const NekDouble>& physvals)
1075  {
1076  return v_PhysEvaluate(coords,physvals);
1077  }
1078 
1079 
1080  /** \brief This function evaluates the expansion at a single
1081  * (arbitrary) point of the domain
1082  *
1083  * This function is a wrapper around the virtual function
1084  * \a v_PhysEvaluate()
1085  *
1086  * Based on the value of the expansion at the quadrature
1087  * points provided in \a physvals, this function
1088  * calculates the value of the expansion at an arbitrary
1089  * single points associated with the interpolation
1090  * matrices provided in \f$ I \f$.
1091  *
1092  * \param I an Array of lagrange interpolantes evaluated
1093  * at the coordinate and going through the local physical
1094  * quadrature
1095  * \param physvals the interpolated field at the quadrature points
1096  *
1097  * \return returns the value of the expansion at the
1098  * single point
1099  */
1100  NekDouble PhysEvaluate(const Array<OneD, DNekMatSharedPtr>& I,
1101  const Array<OneD, const NekDouble >& physvals)
1102  {
1103  return v_PhysEvaluate(I,physvals);
1104  }
1105 
1106 
1107  /**
1108  * \brief Convert local cartesian coordinate \a xi into local
1109  * collapsed coordinates \a eta
1110  **/
1111  void LocCoordToLocCollapsed(const Array<OneD, const NekDouble>& xi,
1112  Array<OneD, NekDouble>& eta)
1113  {
1114  v_LocCoordToLocCollapsed(xi,eta);
1115  }
1116 
1117  const boost::shared_ptr<SpatialDomains::GeomFactors>& GetMetricInfo(void) const
1118  {
1119  return v_GetMetricInfo();
1120  }
1121 
1122  /// \brief Get the element id of this expansion when used
1123  /// in a list by returning value of #m_elmt_id
1124  STD_REGIONS_EXPORT virtual int v_GetElmtId();
1125 
1126  STD_REGIONS_EXPORT virtual const Array<OneD, const NekDouble>& v_GetPhysNormals(void);
1127 
1128  STD_REGIONS_EXPORT virtual void v_SetPhysNormals(Array<OneD, const NekDouble> &normal);
1129 
1130  STD_REGIONS_EXPORT virtual void v_SetUpPhysNormals(const int edge);
1131 
1132  STD_REGIONS_EXPORT virtual int v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset);
1133 
1134  /**
1135  * @brief Unpack data from input file assuming it comes from the
1136  * same expansion type.
1137  * @see StdExpansion::ExtractDataToCoeffs
1138  */
1139  STD_REGIONS_EXPORT virtual void v_ExtractDataToCoeffs(const NekDouble *data,
1140  const std::vector<unsigned int > &nummodes,
1141  const int nmode_offset,
1142  NekDouble *coeffs);
1143 
1144  STD_REGIONS_EXPORT virtual void v_NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, Array< OneD, NekDouble> &outarray);
1145 
1146  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);
1147 
1149 
1151 
1152 
1154 
1156 
1158 
1160 
1161  /** \brief Function to evaluate the discrete \f$ L_\infty\f$
1162  * error \f$ |\epsilon|_\infty = \max |u - u_{exact}|\f$ where \f$
1163  * u_{exact}\f$ is given by the array \a sol.
1164  *
1165  * This function takes the physical value space array \a m_phys as
1166  * approximate solution
1167  *
1168  * \param sol array of solution function at physical quadrature
1169  * points
1170  * \return returns the \f$ L_\infty \f$ error as a NekDouble.
1171  */
1172  STD_REGIONS_EXPORT NekDouble Linf(const Array<OneD, const NekDouble>& phys, const Array<OneD, const NekDouble>& sol = NullNekDouble1DArray);
1173 
1174  /** \brief Function to evaluate the discrete \f$ L_2\f$ error,
1175  * \f$ | \epsilon |_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2
1176  * dx \right]^{1/2} d\xi_1 \f$ where \f$ u_{exact}\f$ is given by
1177  * the array \a sol.
1178  *
1179  * This function takes the physical value space array \a m_phys as
1180  * approximate solution
1181  *
1182  * \param sol array of solution function at physical quadrature
1183  * points
1184  * \return returns the \f$ L_2 \f$ error as a double.
1185  */
1186  STD_REGIONS_EXPORT NekDouble L2(const Array<OneD, const NekDouble>& phys, const Array<OneD, const NekDouble>& sol = NullNekDouble1DArray);
1187 
1188  /** \brief Function to evaluate the discrete \f$ H^1\f$
1189  * error, \f$ | \epsilon |^1_{2} = \left [ \int^1_{-1} [u -
1190  * u_{exact}]^2 + \nabla(u - u_{exact})\cdot\nabla(u -
1191  * u_{exact})\cdot dx \right]^{1/2} d\xi_1 \f$ where \f$
1192  * u_{exact}\f$ is given by the array \a sol.
1193  *
1194  * This function takes the physical value space array
1195  * \a m_phys as approximate solution
1196  *
1197  * \param sol array of solution function at physical quadrature
1198  * points
1199  * \return returns the \f$ H_1 \f$ error as a double.
1200  */
1201  STD_REGIONS_EXPORT NekDouble H1(const Array<OneD, const NekDouble>& phys, const Array<OneD, const NekDouble>& sol = NullNekDouble1DArray);
1202 
1203  // I/O routines
1204  const NormalVector & GetEdgeNormal(const int edge) const
1205  {
1206  return v_GetEdgeNormal(edge);
1207  }
1208 
1209  void ComputeEdgeNormal(const int edge)
1210  {
1211  v_ComputeEdgeNormal(edge);
1212  }
1213 
1214  void NegateEdgeNormal(const int edge)
1215  {
1216  v_NegateEdgeNormal(edge);
1217  }
1218 
1219  bool EdgeNormalNegated(const int edge)
1220  {
1221  return v_EdgeNormalNegated(edge);
1222  }
1223 
1224  void ComputeFaceNormal(const int face)
1225  {
1226  v_ComputeFaceNormal(face);
1227  }
1228 
1229  void NegateFaceNormal(const int face)
1230  {
1231  v_NegateFaceNormal(face);
1232  }
1233 
1234  void ComputeVertexNormal(const int vertex)
1235  {
1236  v_ComputeVertexNormal(vertex);
1237  }
1238 
1239  const NormalVector & GetFaceNormal(const int face) const
1240  {
1241  return v_GetFaceNormal(face);
1242  }
1243 
1244  const NormalVector & GetVertexNormal(const int vertex) const
1245  {
1246  return v_GetVertexNormal(vertex);
1247  }
1248 
1249  const NormalVector & GetSurfaceNormal(const int id) const
1250  {
1251  return v_GetSurfaceNormal(id);
1252  }
1253 
1255  {
1257  for (int i = 0; i < m_base.num_elements(); ++i)
1258  {
1259  p.push_back(m_base[i]->GetPointsKey());
1260  }
1261  return p;
1262  }
1263 
1264  STD_REGIONS_EXPORT Array<OneD, unsigned int>
1266  {
1267  return v_GetEdgeInverseBoundaryMap(eid);
1268  }
1269 
1270  STD_REGIONS_EXPORT Array<OneD, unsigned int>
1272  {
1273  return v_GetFaceInverseBoundaryMap(fid,faceOrient);
1274  }
1275 
1277  const DNekScalMatSharedPtr & m_transformationmatrix)
1278  {
1280  m_transformationmatrix);
1281  }
1282 
1283 
1284  /** \brief This function performs an interpolation from
1285  * the physical space points provided at input into an
1286  * array of equispaced points which are not the collapsed
1287  * coordinate. So for a tetrahedron you will only get a
1288  * tetrahedral number of values.
1289  *
1290  * This is primarily used for output purposes to get a
1291  * better distribution of points more suitable for most
1292  * postprocessing
1293  */
1295  const Array<OneD, const NekDouble> &inarray,
1296  Array<OneD, NekDouble> &outarray);
1297 
1298  /** \brief This function provides the connectivity of
1299  * local simplices (triangles or tets) to connect the
1300  * equispaced data points provided by
1301  * PhysInterpToSimplexEquiSpaced
1302  *
1303  * This is a virtual call to the function
1304  * \a v_GetSimplexEquiSpaceConnectivity
1305  */
1307  Array<OneD, int> &conn,
1308  bool standard = true)
1309  {
1310  v_GetSimplexEquiSpacedConnectivity(conn,standard);
1311  }
1312 
1313  template<class T>
1314  boost::shared_ptr<T> as()
1315  {
1316 #if defined __INTEL_COMPILER && BOOST_VERSION > 105200
1317  typedef typename boost::shared_ptr<T>::element_type E;
1318  E * p = dynamic_cast< E* >( shared_from_this().get() );
1319  ASSERTL1(p, "Cannot perform cast");
1320  return boost::shared_ptr<T>( shared_from_this(), p );
1321 #else
1322  return boost::dynamic_pointer_cast<T>( shared_from_this() );
1323 #endif
1324  }
1325 
1326  protected:
1327  Array<OneD, LibUtilities::BasisSharedPtr> m_base; /**< Bases needed for the expansion */
1329  int m_ncoeffs; /**< Total number of coefficients used in the expansion */
1333 
1335  {
1336  return v_CreateStdMatrix(mkey);
1337  }
1338 
1339  /** \brief Create the static condensation of a matrix when
1340  using a boundary interior decomposition
1341 
1342  If a matrix system can be represented by
1343  \f$ Mat = \left [ \begin{array}{cc}
1344  A & B \\
1345  C & D \end{array} \right ] \f$
1346  This routine creates a matrix containing the statically
1347  condense system of the form
1348  \f$ Mat = \left [ \begin{array}{cc}
1349  A - B D^{-1} C & B D^{-1} \\
1350  D^{-1} C & D^{-1} \end{array} \right ] \f$
1351  **/
1353 
1354  /** \brief Create an IndexMap which contains mapping information linking any specific
1355  element shape with either its boundaries, edges, faces, verteces, etc.
1356 
1357  The index member of the IndexMapValue struct gives back an integer associated with an entity index
1358  The sign member of the same struct gives back a sign to algebrically apply entities orientation
1359  **/
1361 
1362  STD_REGIONS_EXPORT void BwdTrans_MatOp(const Array<OneD, const NekDouble>& inarray,
1363  Array<OneD, NekDouble> &outarray);
1364 
1365  void BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
1366  Array<OneD, NekDouble> &outarray)
1367  {
1368  v_BwdTrans_SumFac(inarray,outarray);
1369  }
1370 
1371  void IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
1372  Array<OneD, NekDouble> &outarray)
1373  {
1374  v_IProductWRTBase_SumFac(inarray,outarray);
1375  }
1376 
1377  void IProductWRTDerivBase_SumFac(const int dir,
1378  const Array<OneD, const NekDouble>& inarray,
1379  Array<OneD, NekDouble> &outarray)
1380  {
1381  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
1382  }
1383 
1384  // The term _MatFree denotes that the action of the MatrixOperation
1385  // is done withouth actually using the matrix (which then needs to be stored/calculated).
1386  // Although this does not strictly mean that no matrix operations are involved in the
1387  // evaluation of the operation, we use this term in the same context used as in the following
1388  // paper:
1389  // R. C. Kirby, M. G. Knepley, A. Logg, and L. R. Scott,
1390  // "Optimizing the evaluation of finite element matrices," SISC 27:741-758 (2005)
1391  STD_REGIONS_EXPORT void GeneralMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1392  Array<OneD,NekDouble> &outarray,
1393  const StdMatrixKey &mkey);
1394 
1395  STD_REGIONS_EXPORT void MassMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1396  Array<OneD,NekDouble> &outarray,
1397  const StdMatrixKey &mkey);
1398 
1399  void LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1400  Array<OneD,NekDouble> &outarray,
1401  const StdMatrixKey &mkey)
1402  {
1403  v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1404  }
1405 
1407  const Array<OneD, const NekDouble> &inarray,
1408  Array<OneD, NekDouble> &outarray,
1409  Array<OneD, NekDouble> &wsp)
1410  {
1411  v_LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp);
1412  }
1413 
1414  STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree_GenericImpl(const Array<OneD, const NekDouble> &inarray,
1415  Array<OneD,NekDouble> &outarray,
1416  const StdMatrixKey &mkey);
1417 
1418  STD_REGIONS_EXPORT void LaplacianMatrixOp_MatFree(const int k1, const int k2,
1419  const Array<OneD, const NekDouble> &inarray,
1420  Array<OneD,NekDouble> &outarray,
1421  const StdMatrixKey &mkey);
1422 
1424  const Array<OneD, const NekDouble> &inarray,
1425  Array<OneD,NekDouble> &outarray,
1426  const StdMatrixKey &mkey);
1427 
1428  STD_REGIONS_EXPORT void WeakDirectionalDerivMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1429  Array<OneD,NekDouble> &outarray,
1430  const StdMatrixKey &mkey);
1431 
1432  STD_REGIONS_EXPORT void MassLevelCurvatureMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1433  Array<OneD,NekDouble> &outarray,
1434  const StdMatrixKey &mkey);
1435 
1436  STD_REGIONS_EXPORT void LinearAdvectionDiffusionReactionMatrixOp_MatFree( const Array<OneD, const NekDouble> &inarray,
1437  Array<OneD,NekDouble> &outarray,
1438  const StdMatrixKey &mkey,
1439  bool addDiffusionTerm = true);
1440 
1441  void HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1442  Array<OneD,NekDouble> &outarray,
1443  const StdMatrixKey &mkey)
1444  {
1445  v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1446  }
1447 
1448  STD_REGIONS_EXPORT void HelmholtzMatrixOp_MatFree_GenericImpl(const Array<OneD, const NekDouble> &inarray,
1449  Array<OneD,NekDouble> &outarray,
1450  const StdMatrixKey &mkey);
1451 
1453  Array<OneD, const NekDouble> &inarray,
1454  Array<OneD, NekDouble> &outarray);
1455 
1457  Array<OneD, NekDouble> &coeffs,
1459 
1461  const Array<OneD, const NekDouble> &Lcoord,
1462  const Array<OneD, const NekDouble> &physvals);
1463 
1464  private:
1465  // Virtual functions
1466  STD_REGIONS_EXPORT virtual int v_GetNverts() const = 0;
1467  STD_REGIONS_EXPORT virtual int v_GetNedges() const;
1468 
1469  STD_REGIONS_EXPORT virtual int v_GetNfaces() const;
1470 
1471  STD_REGIONS_EXPORT virtual int v_NumBndryCoeffs() const;
1472 
1473  STD_REGIONS_EXPORT virtual int v_NumDGBndryCoeffs() const;
1474 
1475  STD_REGIONS_EXPORT virtual int v_GetEdgeNcoeffs(const int i) const;
1476 
1477  STD_REGIONS_EXPORT virtual int v_GetTotalEdgeIntNcoeffs() const;
1478 
1479  STD_REGIONS_EXPORT virtual int v_GetEdgeNumPoints(const int i) const;
1480 
1481  STD_REGIONS_EXPORT virtual int v_DetCartesianDirOfEdge(const int edge);
1482 
1483  STD_REGIONS_EXPORT virtual const LibUtilities::BasisKey v_DetEdgeBasisKey(const int i) const;
1484 
1485  STD_REGIONS_EXPORT virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const;
1486 
1487  STD_REGIONS_EXPORT virtual int v_GetFaceNumPoints(const int i) const;
1488 
1489  STD_REGIONS_EXPORT virtual int v_GetFaceNcoeffs(const int i) const;
1490 
1491  STD_REGIONS_EXPORT virtual int v_GetFaceIntNcoeffs(const int i) const;
1492 
1493  STD_REGIONS_EXPORT virtual int v_GetTotalFaceIntNcoeffs() const;
1494 
1495  STD_REGIONS_EXPORT virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const;
1496 
1498 
1500 
1501  STD_REGIONS_EXPORT virtual int v_GetShapeDimension() const;
1502 
1504 
1505  STD_REGIONS_EXPORT virtual void v_BwdTrans (const Array<OneD, const NekDouble>& inarray,
1506  Array<OneD, NekDouble> &outarray) = 0;
1507 
1508  /**
1509  * @brief Transform a given function from physical quadrature space
1510  * to coefficient space.
1511  * @see StdExpansion::FwdTrans
1512  */
1513  STD_REGIONS_EXPORT virtual void v_FwdTrans (
1514  const Array<OneD, const NekDouble>& inarray,
1515  Array<OneD, NekDouble> &outarray) = 0;
1516 
1517  /**
1518  * @brief Calculates the inner product of a given function \a f
1519  * with the different modes of the expansion
1520  */
1522  const Array<OneD, const NekDouble>& inarray,
1523  Array<OneD, NekDouble> &outarray) = 0;
1524 
1526  const Array<OneD, const NekDouble>& base,
1527  const Array<OneD, const NekDouble>& inarray,
1528  Array<OneD, NekDouble>& outarray,
1529  int coll_check)
1530  {
1531  ASSERTL0(false, "StdExpansion::v_IProductWRTBase has no (and should have no) implementation");
1532  }
1533 
1534  STD_REGIONS_EXPORT virtual void v_IProductWRTDerivBase (const int dir,
1535  const Array<OneD, const NekDouble>& inarray,
1536  Array<OneD, NekDouble> &outarray);
1537 
1538  STD_REGIONS_EXPORT virtual void v_FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
1539  Array<OneD, NekDouble> &outarray);
1540 
1541  STD_REGIONS_EXPORT virtual NekDouble v_Integral(const Array<OneD, const NekDouble>& inarray );
1542 
1543  STD_REGIONS_EXPORT virtual void v_PhysDeriv (const Array<OneD, const NekDouble>& inarray,
1544  Array<OneD, NekDouble> &out_d1,
1545  Array<OneD, NekDouble> &out_d2,
1546  Array<OneD, NekDouble> &out_d3);
1547 
1548  STD_REGIONS_EXPORT virtual void v_PhysDeriv_s (const Array<OneD, const NekDouble>& inarray,
1549  Array<OneD, NekDouble> &out_ds);
1550 
1551  STD_REGIONS_EXPORT virtual void v_PhysDeriv_n(const Array<OneD, const NekDouble>& inarray,
1552  Array<OneD, NekDouble>& out_dn);
1553  STD_REGIONS_EXPORT virtual void v_PhysDeriv(const int dir,
1554  const Array<OneD, const NekDouble>& inarray,
1555  Array<OneD, NekDouble> &out_d0);
1556 
1557  STD_REGIONS_EXPORT virtual void v_PhysDirectionalDeriv(const Array<OneD, const NekDouble>& inarray,
1558  const Array<OneD, const NekDouble>& direction,
1559  Array<OneD, NekDouble> &outarray);
1560 
1561  STD_REGIONS_EXPORT virtual void v_StdPhysDeriv (const Array<OneD, const NekDouble>& inarray,
1562  Array<OneD, NekDouble> &out_d1,
1563  Array<OneD, NekDouble> &out_d2,
1564  Array<OneD, NekDouble> &out_d3);
1565 
1566  STD_REGIONS_EXPORT virtual void v_StdPhysDeriv (const int dir,
1567  const Array<OneD, const NekDouble>& inarray,
1568  Array<OneD, NekDouble> &outarray);
1569 
1570  STD_REGIONS_EXPORT virtual void v_AddRobinMassMatrix(const int edgeid, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat);
1571 
1572  STD_REGIONS_EXPORT virtual void v_AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble> &primCoeffs, Array<OneD, NekDouble> &coeffs);
1573 
1574  STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
1575 
1576  STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, DNekMatSharedPtr >& I, const Array<OneD, const NekDouble> & physvals);
1577 
1579  const Array<OneD, const NekDouble>& xi,
1580  Array<OneD, NekDouble>& eta);
1581 
1582 
1583  STD_REGIONS_EXPORT virtual void v_FillMode(const int mode, Array<OneD, NekDouble> &outarray);
1584 
1586 
1588 
1589  STD_REGIONS_EXPORT virtual void v_GetCoords(Array<OneD, NekDouble> &coords_0,
1590  Array<OneD, NekDouble> &coords_1,
1591  Array<OneD, NekDouble> &coords_2);
1592 
1593  STD_REGIONS_EXPORT virtual void v_GetCoord(const Array<OneD, const NekDouble>& Lcoord,
1594  Array<OneD, NekDouble> &coord);
1595 
1596  STD_REGIONS_EXPORT virtual int v_GetCoordim(void);
1597 
1598  STD_REGIONS_EXPORT virtual void v_GetBoundaryMap(Array<OneD, unsigned int>& outarray);
1599 
1600  STD_REGIONS_EXPORT virtual void v_GetInteriorMap(Array<OneD, unsigned int>& outarray);
1601 
1602  STD_REGIONS_EXPORT virtual int v_GetVertexMap(int localVertexId,
1603  bool useCoeffPacking = false);
1604 
1605  STD_REGIONS_EXPORT virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient,
1606  Array<OneD, unsigned int> &maparray,
1607  Array<OneD, int> &signarray);
1608 
1609  STD_REGIONS_EXPORT virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient,
1610  Array<OneD, unsigned int> &maparray,
1611  Array<OneD, int> &signarray);
1612 
1613  STD_REGIONS_EXPORT virtual void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient,
1614  Array<OneD, unsigned int> &maparray,
1615  Array<OneD, int> &signarray);
1616 
1617  STD_REGIONS_EXPORT virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient,
1618  Array<OneD, unsigned int> &maparray,
1619  Array<OneD, int> &signarray,
1620  int nummodesA = -1, int nummodesB = -1);
1621 
1622  /**
1623  * @brief Extract the physical values along edge \a edge from \a
1624  * inarray into \a outarray following the local edge orientation
1625  * and point distribution defined by defined in \a EdgeExp.
1626  */
1627  STD_REGIONS_EXPORT virtual void v_GetEdgePhysVals(const int edge, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray);
1628 
1629  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);
1630 
1631  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);
1632 
1633  STD_REGIONS_EXPORT virtual void v_GetVertexPhysVals(const int vertex, const Array<OneD, const NekDouble> &inarray, NekDouble &outarray);
1634 
1635  STD_REGIONS_EXPORT virtual void v_GetEdgeInterpVals(const int edge,
1636  const Array<OneD, const NekDouble> &inarray,Array<OneD,NekDouble> &outarray);
1637 
1639  const int edge,
1640  Array<OneD, NekDouble> &outarray);
1641 
1643  const int face,
1644  const boost::shared_ptr<StdExpansion> &FaceExp,
1645  const Array<OneD, const NekDouble> &inarray,
1646  Array<OneD, NekDouble> &outarray,
1647  StdRegions::Orientation orient);
1648 
1650  const Array<OneD, const NekDouble> &inarray,
1651  Array<OneD, NekDouble> &outarray);
1652 
1654  const Array<OneD, const NekDouble> &inarray,
1655  Array<OneD, NekDouble> &outarray);
1656 
1657  STD_REGIONS_EXPORT virtual const boost::shared_ptr<SpatialDomains::GeomFactors>& v_GetMetricInfo() const;
1658 
1659  STD_REGIONS_EXPORT virtual void v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
1660  Array<OneD, NekDouble> &outarray);
1661 
1662  STD_REGIONS_EXPORT virtual void v_IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
1663  Array<OneD, NekDouble> &outarray);
1664 
1665  STD_REGIONS_EXPORT virtual void v_IProductWRTDerivBase_SumFac(const int dir,
1666  const Array<OneD, const NekDouble>& inarray,
1667  Array<OneD, NekDouble> &outarray);
1668 
1669  STD_REGIONS_EXPORT virtual void v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
1670  Array<OneD,NekDouble> &outarray,
1671  const StdMatrixKey &mkey);
1672 
1673  STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
1674  Array<OneD,NekDouble> &outarray,
1675  const StdMatrixKey &mkey);
1676 
1677  STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD,NekDouble> &array,
1678  const StdMatrixKey &mkey);
1679 
1681  int numMin,
1682  const Array<OneD, const NekDouble> &inarray,
1683  Array<OneD,NekDouble> &outarray);
1684 
1685  STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp(const int k1, const int k2,
1686  const Array<OneD, const NekDouble> &inarray,
1687  Array<OneD,NekDouble> &outarray,
1688  const StdMatrixKey &mkey);
1689 
1690  STD_REGIONS_EXPORT virtual void v_WeakDerivMatrixOp(const int i,
1691  const Array<OneD, const NekDouble> &inarray,
1692  Array<OneD,NekDouble> &outarray,
1693  const StdMatrixKey &mkey);
1694 
1695  STD_REGIONS_EXPORT virtual void v_WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
1696  Array<OneD,NekDouble> &outarray,
1697  const StdMatrixKey &mkey);
1698 
1699  STD_REGIONS_EXPORT virtual void v_MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
1700  Array<OneD,NekDouble> &outarray,
1701  const StdMatrixKey &mkey);
1702 
1704  const NekDouble> &inarray,
1705  Array<OneD,NekDouble> &outarray,
1706  const StdMatrixKey &mkey,
1707  bool addDiffusionTerm=true);
1708 
1709  STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
1710  Array<OneD,NekDouble> &outarray,
1711  const StdMatrixKey &mkey);
1712 
1713  STD_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1714  Array<OneD,NekDouble> &outarray,
1715  const StdMatrixKey &mkey);
1716 
1718  const Array<OneD, const NekDouble> &inarray,
1719  Array<OneD, NekDouble> &outarray,
1720  Array<OneD, NekDouble> &wsp);
1721 
1722  STD_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1723  Array<OneD,NekDouble> &outarray,
1724  const StdMatrixKey &mkey);
1725 
1726  STD_REGIONS_EXPORT virtual const NormalVector & v_GetEdgeNormal(const int edge) const;
1727 
1728  STD_REGIONS_EXPORT virtual void v_ComputeEdgeNormal(const int edge);
1729 
1730  STD_REGIONS_EXPORT virtual void v_NegateEdgeNormal(const int edge);
1731 
1732  STD_REGIONS_EXPORT virtual bool v_EdgeNormalNegated(const int edge);
1733 
1734  STD_REGIONS_EXPORT virtual void v_ComputeFaceNormal(const int face);
1735 
1736  STD_REGIONS_EXPORT virtual void v_NegateFaceNormal(const int face);
1737 
1738  STD_REGIONS_EXPORT virtual const NormalVector & v_GetVertexNormal(const int vertex) const;
1739 
1740  STD_REGIONS_EXPORT virtual void v_ComputeVertexNormal(const int vertex);
1741 
1742  STD_REGIONS_EXPORT virtual const NormalVector & v_GetFaceNormal(const int face) const;
1743  STD_REGIONS_EXPORT virtual const NormalVector &
1744  v_GetSurfaceNormal(const int id) const;
1745 
1746  STD_REGIONS_EXPORT virtual Array<OneD, unsigned int>
1747  v_GetEdgeInverseBoundaryMap(int eid);
1748 
1749  STD_REGIONS_EXPORT virtual Array<OneD, unsigned int>
1751 
1753 
1755  Array<OneD, int> &conn,
1756  bool standard = true);
1757  };
1758 
1759 
1760  typedef boost::shared_ptr<StdExpansion> StdExpansionSharedPtr;
1761  typedef std::vector< StdExpansionSharedPtr > StdExpansionVector;
1763 
1764  /**
1765  * This function is a wrapper around the virtual function
1766  * \a v_FwdTrans()
1767  *
1768  * Given a function evaluated at the quadrature points, this
1769  * function calculates the expansion coefficients such that the
1770  * resulting expansion approximates the original function.
1771  *
1772  * The calculation of the expansion coefficients is done using a
1773  * Galerkin projection. This is equivalent to the operation:
1774  * \f[ \mathbf{\hat{u}} = \mathbf{M}^{-1} \mathbf{I}\f]
1775  * where
1776  * - \f$\mathbf{M}[p][q]= \int\phi_p(\mathbf{\xi})\phi_q(
1777  * \mathbf{\xi}) d\mathbf{\xi}\f$ is the Mass matrix
1778  * - \f$\mathbf{I}[p] = \int\phi_p(\mathbf{\xi}) u(\mathbf{\xi})
1779  * d\mathbf{\xi}\f$
1780  *
1781  * This function takes the array \a inarray as the values of the
1782  * function evaluated at the quadrature points
1783  * (i.e. \f$\mathbf{u}\f$),
1784  * and stores the resulting coefficients \f$\mathbf{\hat{u}}\f$
1785  * in the \a outarray
1786  *
1787  * @param inarray array of the function discretely evaluated at the
1788  * quadrature points
1789  *
1790  * @param outarray array of the function coefficieints
1791  */
1792  inline void StdExpansion::FwdTrans (const Array<OneD, const NekDouble>& inarray,
1793  Array<OneD, NekDouble> &outarray)
1794  {
1795  v_FwdTrans(inarray,outarray);
1796  }
1797 
1798  } //end of namespace
1799 } //end of namespace
1800 
1801 #endif //STANDARDDEXPANSION_H