Nektar++
StdNodalPrismExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdNodalPrismExp.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Nodal prismatic routines built upon StdExpansion3D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 
39 namespace Nektar
40 {
41  namespace StdRegions
42  {
44  const LibUtilities::BasisKey &Ba,
45  const LibUtilities::BasisKey &Bb,
46  const LibUtilities::BasisKey &Bc,
48  StdExpansion (LibUtilities::StdPrismData::getNumberOfCoefficients(
49  Ba.GetNumModes(),
50  Bb.GetNumModes(),
51  Bc.GetNumModes()),
52  3,Ba,Bb,Bc),
53  StdExpansion3D(LibUtilities::StdPrismData::getNumberOfCoefficients(
54  Ba.GetNumModes(),
55  Bb.GetNumModes(),
56  Bc.GetNumModes()),
57  Ba,Bb,Bc),
58  StdPrismExp (Ba,Bb,Bc),
59  m_nodalPointsKey(Ba.GetNumModes(),Ntype)
60  {
61  ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
62  "order in 'a' direction is higher than order "
63  "in 'c' direction");
64  }
65 
67  StdExpansion(T),
68  StdExpansion3D(T),
69  StdPrismExp(T),
70  m_nodalPointsKey(T.m_nodalPointsKey)
71  {
72  }
73 
75  {
76  }
77 
78 
79 
81  {
82  return true;
83  }
84 
85  //-------------------------------
86  // Nodal basis specific routines
87  //-------------------------------
88 
90  const Array<OneD, const NekDouble>& inarray,
91  Array<OneD, NekDouble>& outarray)
92  {
96  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
97 
99  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
100  modal = (*inv_vdm) * nodal;
101  }
102 
103  // Operate with transpose of NodalToModal transformation
105  const Array<OneD, const NekDouble>& inarray,
106  Array<OneD, NekDouble>& outarray)
107  {
111  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
112 
113  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eCopy);
114  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
115  modal = Transpose(*inv_vdm) * nodal;
116  }
117 
119  const Array<OneD, const NekDouble>& inarray,
120  Array<OneD, NekDouble>& outarray)
121  {
122  StdMatrixKey Nkey(eNBasisTrans, DetShapeType(), *this,
125  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
126 
127  // Multiply out matrix
128  NekVector<NekDouble> modal(m_ncoeffs,inarray,eWrapper);
129  NekVector<NekDouble> nodal(m_ncoeffs,outarray,eWrapper);
130  nodal = (*vdm)*modal;
131  }
132 
137  {
138  // Get 3D nodal distribution.
139  LibUtilities::PointsManager()[m_nodalPointsKey]->GetPoints(x,y,z);
140  }
141 
143  {
144  int i,j;
147  DNekMatSharedPtr Mat;
148 
151  GetNodalPoints(r,s,t);
152 
153  //Store the values of m_phys in a temporary array
154  int nqtot = GetTotPoints();
155  Array<OneD,NekDouble> tmp_phys(nqtot);
156 
157  for(i = 0; i < m_ncoeffs; ++i)
158  {
159  // fill physical space with mode i
160  StdPrismExp::v_FillMode(i,tmp_phys);
161 
162  // interpolate mode i to the Nodal points 'j' and
163  // store in outarray
164  for(j = 0; j < m_ncoeffs; ++j)
165  {
166  c[0] = r[j];
167  c[1] = s[j];
168  c[2] = t[j];
169  (*Mat)(j,i) = StdPrismExp::v_PhysEvaluate(c,tmp_phys);
170  }
171  }
172 
173  return Mat;
174  }
175 
176 
177  //---------------------------------------
178  // Transforms
179  //---------------------------------------
180 
182  const Array<OneD, const NekDouble>& inarray,
183  Array<OneD, NekDouble>& outarray)
184  {
185  v_BwdTrans_SumFac(inarray,outarray);
186  }
187 
189  const Array<OneD, const NekDouble>& inarray,
190  Array<OneD, NekDouble>& outarray)
191  {
193  NodalToModal(inarray,tmp);
194  StdPrismExp::v_BwdTrans_SumFac(tmp,outarray);
195  }
196 
198  const Array<OneD, const NekDouble>& inarray,
199  Array<OneD, NekDouble>& outarray)
200  {
201  v_IProductWRTBase(inarray,outarray);
202 
203  // get Mass matrix inverse
204  StdMatrixKey masskey(eInvMass, DetShapeType(), *this,
207  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
208 
209  // copy inarray in case inarray == outarray
210  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
212 
213  out = (*matsys)*in;
214  }
215 
216 
217  //---------------------------------------
218  // Inner product functions
219  //---------------------------------------
220 
222  const Array<OneD, const NekDouble>& inarray,
223  Array<OneD, NekDouble>& outarray)
224  {
225  v_IProductWRTBase_SumFac(inarray,outarray);
226  }
227 
229  const Array<OneD, const NekDouble>& inarray,
230  Array<OneD, NekDouble>& outarray,
231  bool multiplybyweights)
232  {
233  StdPrismExp::v_IProductWRTBase_SumFac(inarray,outarray,multiplybyweights);
234  NodalToModalTranspose(outarray,outarray);
235  }
236 
238  const int dir,
239  const Array<OneD, const NekDouble>& inarray,
240  Array<OneD, NekDouble>& outarray)
241  {
242  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
243  }
244 
246  const int dir,
247  const Array<OneD, const NekDouble>& inarray,
248  Array<OneD, NekDouble>& outarray)
249  {
250  StdPrismExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
251  NodalToModalTranspose(outarray,outarray);
252  }
253 
254  //---------------------------------------
255  // Evaluation functions
256  //---------------------------------------
257 
259  const int mode,
260  Array<OneD, NekDouble> &outarray)
261  {
262  ASSERTL2(mode >= m_ncoeffs,
263  "calling argument mode is larger than total expansion order");
264 
265  Vmath::Zero(m_ncoeffs, outarray, 1);
266  outarray[mode] = 1.0;
267  v_BwdTrans(outarray,outarray);
268  }
269 
270 
271  //---------------------------------------
272  // Mapping functions
273  //---------------------------------------
274 
275  /*
276  void StdNodalTriExp::v_GetFaceToElementMap(
277  const int fid,
278  const FaceOrientation faceOrient,
279  Array<OneD, unsigned int> &maparray,
280  Array<OneD, int> &signarray,
281  int nummodesA,
282  int nummodesB)
283  {
284  int P, Q, i, j, k, idx = 0, nFaceCoeffs = 0;
285 
286  ASSERTL0(fid >= 0 && fid <= 3,
287  "Local face ID must be between 0 and 3");
288 
289  if (nummodesA == -1)
290  {
291  switch(fid)
292  {
293  case 0:
294  nummodesA = m_base[0]->GetNumModes();
295  nummodesB = m_base[1]->GetNumModes();
296  break;
297  case 1:
298  nummodesA = m_base[0]->GetNumModes();
299  nummodesB = m_base[2]->GetNumModes();
300  break;
301  case 2:
302  case 3:
303  nummodesA = m_base[1]->GetNumModes();
304  nummodesB = m_base[2]->GetNumModes();
305  break;
306  }
307  }
308 
309  P = nummodesA;
310  Q = nummodesB;
311  nFaceCoeffs = Q + ((P-1)*(1 + 2*(Q-1) - (P-1)))/2;
312 
313  if (maparray.num_elements() != nFaceCoeffs)
314  {
315  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
316  }
317 
318  if (signarray.num_elements() != nFaceCoeffs)
319  {
320  signarray = Array<OneD, int>(nFaceCoeffs,1);
321  }
322  else
323  {
324  fill(signarray.get(), signarray.get()+nFaceCoeffs, 1);
325  }
326 
327  switch(fid)
328  {
329  case 0:
330  // Add vertices.
331  maparray[idx++] = 0;
332  maparray[idx++] = 1;
333  maparray[idx++] = 2;
334 
335  // Add edges.
336  for (i = 2; i < P; ++i)
337  {
338  maparray[idx++] = ;
339  }
340  }
341  }
342  */
343 
344  int StdNodalPrismExp::v_GetVertexMap(const int localVertexId,
345  bool useCoeffPacking)
346  {
347  ASSERTL0(false,"Needs setting up");
348  return localVertexId;
349  }
350 
352  Array<OneD, unsigned int>& outarray)
353  {
354  unsigned int i;
355  const unsigned int nBndryCoeff = NumBndryCoeffs();
356 
357  if (outarray.num_elements() != nBndryCoeff)
358  {
359  outarray = Array<OneD, unsigned int>(nBndryCoeff);
360  }
361 
362  for (i = 0; i < nBndryCoeff; i++)
363  {
364  outarray[i] = i;
365  }
366  }
367 
369  Array<OneD, unsigned int>& outarray)
370  {
371  unsigned int i;
372  const unsigned int nBndryCoeff = NumBndryCoeffs();
373 
374  if (outarray.num_elements() != m_ncoeffs-nBndryCoeff)
375  {
376  outarray = Array<OneD, unsigned int>(
377  m_ncoeffs-nBndryCoeff);
378  }
379 
380  for (i = nBndryCoeff; i < m_ncoeffs; i++)
381  {
382  outarray[i-nBndryCoeff] = i;
383  }
384  }
385 
386 
387  //---------------------------------------
388  // Wrapper functions
389  //---------------------------------------
390 
392  {
393  DNekMatSharedPtr Mat;
394 
395  switch(mkey.GetMatrixType())
396  {
397  case eNBasisTrans:
398  Mat = GenNBasisTransMatrix();
399  break;
400  default:
402  break;
403  }
404 
405  return Mat;
406  }
407 
409  const StdMatrixKey &mkey)
410  {
411  return StdNodalPrismExp::v_GenMatrix(mkey);
412  }
413  } // end of namespace
414 } // end of namespace
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Inner product of inarray over region with respect to the object's default expansion basis; output in ...
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:454
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
void GetNodalPoints(Array< OneD, const NekDouble > &x, Array< OneD, const NekDouble > &y, Array< OneD, const NekDouble > &z)
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool mult=true)
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
The base class for all shapes.
Definition: StdExpansion.h:69
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void NodalToModal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Class representing a prismatic element in reference space.
Definition: StdPrismExp.h:49
PointsManagerT & PointsManager(void)
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the basis B=base0*base1*base2 and put into out...
LibUtilities::PointsKey m_nodalPointsKey
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
virtual int v_GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Describes the specification for a Basis.
Definition: Basis.h:50
PointsType GetPointsType() const
Definition: Points.h:111
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:226
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:249