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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Nodal prismatic routines built upon StdExpansion3D
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
39 
40 namespace Nektar
41 {
42 namespace StdRegions
43 {
45  const LibUtilities::BasisKey &Bb,
46  const LibUtilities::BasisKey &Bc,
48  : StdExpansion(LibUtilities::StdPrismData::getNumberOfCoefficients(
49  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
50  3, Ba, Bb, Bc),
51  StdExpansion3D(LibUtilities::StdPrismData::getNumberOfCoefficients(
52  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
53  Ba, Bb, Bc),
54  StdPrismExp(Ba, Bb, Bc), m_nodalPointsKey(Ba.GetNumModes(), Ntype)
55 {
56  ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
57  "order in 'a' direction is higher than order "
58  "in 'c' direction");
59 }
60 
63  m_nodalPointsKey(T.m_nodalPointsKey)
64 {
65 }
66 
68 {
69 }
70 
72 {
73  return true;
74 }
75 
76 //-------------------------------
77 // Nodal basis specific routines
78 //-------------------------------
79 
81  Array<OneD, NekDouble> &outarray)
82 {
86  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
87 
88  NekVector<NekDouble> nodal(m_ncoeffs, inarray, eWrapper);
89  NekVector<NekDouble> modal(m_ncoeffs, outarray, eWrapper);
90  modal = (*inv_vdm) * nodal;
91 }
92 
93 // Operate with transpose of NodalToModal transformation
95  const Array<OneD, const NekDouble> &inarray,
96  Array<OneD, NekDouble> &outarray)
97 {
101  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
102 
103  NekVector<NekDouble> nodal(m_ncoeffs, inarray, eCopy);
104  NekVector<NekDouble> modal(m_ncoeffs, outarray, eWrapper);
105  modal = Transpose(*inv_vdm) * nodal;
106 }
107 
109  Array<OneD, NekDouble> &outarray)
110 {
113  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
114 
115  // Multiply out matrix
116  NekVector<NekDouble> modal(m_ncoeffs, inarray, eWrapper);
117  NekVector<NekDouble> nodal(m_ncoeffs, outarray, eWrapper);
118  nodal = (*vdm) * modal;
119 }
120 
124 {
125  // Get 3D nodal distribution.
126  LibUtilities::PointsManager()[m_nodalPointsKey]->GetPoints(x, y, z);
127 }
128 
130 {
131  int i, j;
134  DNekMatSharedPtr Mat;
135 
137  GetNodalPoints(r, s, t);
138 
139  // Store the values of m_phys in a temporary array
140  int nqtot = GetTotPoints();
141  Array<OneD, NekDouble> tmp_phys(nqtot);
142 
143  for (i = 0; i < m_ncoeffs; ++i)
144  {
145  // fill physical space with mode i
146  StdPrismExp::v_FillMode(i, tmp_phys);
147 
148  // interpolate mode i to the Nodal points 'j' and
149  // store in outarray
150  for (j = 0; j < m_ncoeffs; ++j)
151  {
152  c[0] = r[j];
153  c[1] = s[j];
154  c[2] = t[j];
155  (*Mat)(j, i) = StdPrismExp::v_PhysEvaluate(c, tmp_phys);
156  }
157  }
158 
159  return Mat;
160 }
161 
162 //---------------------------------------
163 // Transforms
164 //---------------------------------------
165 
167  Array<OneD, NekDouble> &outarray)
168 {
169  v_BwdTrans_SumFac(inarray, outarray);
170 }
171 
173  const Array<OneD, const NekDouble> &inarray,
174  Array<OneD, NekDouble> &outarray)
175 {
177  NodalToModal(inarray, tmp);
178  StdPrismExp::v_BwdTrans_SumFac(tmp, outarray);
179 }
180 
182  Array<OneD, NekDouble> &outarray)
183 {
184  v_IProductWRTBase(inarray, outarray);
185 
186  // get Mass matrix inverse
189  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
190 
191  // copy inarray in case inarray == outarray
192  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
193  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
194 
195  out = (*matsys) * in;
196 }
197 
198 //---------------------------------------
199 // Inner product functions
200 //---------------------------------------
201 
203  const Array<OneD, const NekDouble> &inarray,
204  Array<OneD, NekDouble> &outarray)
205 {
206  v_IProductWRTBase_SumFac(inarray, outarray);
207 }
208 
210  const Array<OneD, const NekDouble> &inarray,
211  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
212 {
213  StdPrismExp::v_IProductWRTBase_SumFac(inarray, outarray, multiplybyweights);
214  NodalToModalTranspose(outarray, outarray);
215 }
216 
218  const int dir, const Array<OneD, const NekDouble> &inarray,
219  Array<OneD, NekDouble> &outarray)
220 {
221  v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
222 }
223 
225  const int dir, const Array<OneD, const NekDouble> &inarray,
226  Array<OneD, NekDouble> &outarray)
227 {
228  StdPrismExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
229  NodalToModalTranspose(outarray, outarray);
230 }
231 
232 //---------------------------------------
233 // Evaluation functions
234 //---------------------------------------
235 
236 void StdNodalPrismExp::v_FillMode(const int mode,
237  Array<OneD, NekDouble> &outarray)
238 {
239  ASSERTL2(mode >= m_ncoeffs,
240  "calling argument mode is larger than total expansion order");
241 
242  Vmath::Zero(m_ncoeffs, outarray, 1);
243  outarray[mode] = 1.0;
244  v_BwdTrans(outarray, outarray);
245 }
246 
247 //---------------------------------------
248 // Mapping functions
249 //---------------------------------------
250 
251 /*
252 void StdNodalTriExp::v_GetFaceToElementMap(
253  const int fid,
254  const FaceOrientation faceOrient,
255  Array<OneD, unsigned int> &maparray,
256  Array<OneD, int> &signarray,
257  int nummodesA,
258  int nummodesB)
259 {
260  int P, Q, i, j, k, idx = 0, nFaceCoeffs = 0;
261 
262  ASSERTL0(fid >= 0 && fid <= 3,
263  "Local face ID must be between 0 and 3");
264 
265  if (nummodesA == -1)
266  {
267  switch(fid)
268  {
269  case 0:
270  nummodesA = m_base[0]->GetNumModes();
271  nummodesB = m_base[1]->GetNumModes();
272  break;
273  case 1:
274  nummodesA = m_base[0]->GetNumModes();
275  nummodesB = m_base[2]->GetNumModes();
276  break;
277  case 2:
278  case 3:
279  nummodesA = m_base[1]->GetNumModes();
280  nummodesB = m_base[2]->GetNumModes();
281  break;
282  }
283  }
284 
285  P = nummodesA;
286  Q = nummodesB;
287  nFaceCoeffs = Q + ((P-1)*(1 + 2*(Q-1) - (P-1)))/2;
288 
289  if (maparray.size() != nFaceCoeffs)
290  {
291  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
292  }
293 
294  if (signarray.size() != nFaceCoeffs)
295  {
296  signarray = Array<OneD, int>(nFaceCoeffs,1);
297  }
298  else
299  {
300  fill(signarray.get(), signarray.get()+nFaceCoeffs, 1);
301  }
302 
303  switch(fid)
304  {
305  case 0:
306  // Add vertices.
307  maparray[idx++] = 0;
308  maparray[idx++] = 1;
309  maparray[idx++] = 2;
310 
311  // Add edges.
312  for (i = 2; i < P; ++i)
313  {
314  maparray[idx++] = ;
315  }
316  }
317 }
318 */
319 
320 int StdNodalPrismExp::v_GetVertexMap(const int localVertexId,
321  bool useCoeffPacking)
322 {
323  boost::ignore_unused(useCoeffPacking);
324  ASSERTL0(false, "Needs setting up");
325  return localVertexId;
326 }
327 
329 {
330  unsigned int i;
331  const unsigned int nBndryCoeff = NumBndryCoeffs();
332 
333  if (outarray.size() != nBndryCoeff)
334  {
335  outarray = Array<OneD, unsigned int>(nBndryCoeff);
336  }
337 
338  for (i = 0; i < nBndryCoeff; i++)
339  {
340  outarray[i] = i;
341  }
342 }
343 
345 {
346  unsigned int i;
347  const unsigned int nBndryCoeff = NumBndryCoeffs();
348 
349  if (outarray.size() != m_ncoeffs - nBndryCoeff)
350  {
351  outarray = Array<OneD, unsigned int>(m_ncoeffs - nBndryCoeff);
352  }
353 
354  for (i = nBndryCoeff; i < m_ncoeffs; i++)
355  {
356  outarray[i - nBndryCoeff] = i;
357  }
358 }
359 
360 //---------------------------------------
361 // Wrapper functions
362 //---------------------------------------
363 
365 {
366  DNekMatSharedPtr Mat;
367 
368  switch (mkey.GetMatrixType())
369  {
370  case eNBasisTrans:
371  Mat = GenNBasisTransMatrix();
372  break;
373  default:
375  break;
376  }
377 
378  return Mat;
379 }
380 
382 {
383  return StdNodalPrismExp::v_GenMatrix(mkey);
384 }
385 } // namespace StdRegions
386 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
Describes the specification for a Basis.
Definition: Basis.h:50
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:83
PointsType GetPointsType() const
Definition: Points.h:109
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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.
The base class for all shapes.
Definition: StdExpansion.h:71
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
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...
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...
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::PointsKey m_nodalPointsKey
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
virtual int v_GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
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 ...
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool mult=true)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
void GetNodalPoints(Array< OneD, const NekDouble > &x, Array< OneD, const NekDouble > &y, Array< OneD, const NekDouble > &z)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
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)
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Class representing a prismatic element in reference space.
Definition: StdPrismExp.h:49
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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)
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
PointsManagerT & PointsManager(void)
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:283
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:241
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492