Nektar++
StdNodalTetExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdNodalTetExp.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 tetrahedral 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::StdTetData::getNumberOfCoefficients(
49  Ba.GetNumModes(),
50  Bb.GetNumModes(),
51  Bc.GetNumModes()),
52  3,Ba,Bb,Bc),
53  StdExpansion3D(LibUtilities::StdTetData::getNumberOfCoefficients(
54  Ba.GetNumModes(),
55  Bb.GetNumModes(),
56  Bc.GetNumModes()),
57  Ba,Bb,Bc),
58  StdTetExp (Ba,Bb,Bc),
59  m_nodalPointsKey(Ba.GetNumModes(),Ntype)
60  {
61  ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
62  "order in 'a' direction is higher than order "
63  "in 'b' direction");
64  ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
65  "order in 'a' direction is higher than order "
66  "in 'c' direction");
67  ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(),
68  "order in 'b' direction is higher than order "
69  "in 'c' direction");
70  }
71 
73  StdExpansion(T),
74  StdExpansion3D(T),
75  StdTetExp(T),
76  m_nodalPointsKey(T.m_nodalPointsKey)
77  {
78  }
79 
81  {
82  }
83 
85  {
86  return true;
87  }
88 
89  //-------------------------------
90  // Nodal basis specific routines
91  //-------------------------------
92 
94  const Array<OneD, const NekDouble>& inarray,
95  Array<OneD, NekDouble>& outarray)
96  {
100  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
101 
102  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eWrapper);
103  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
104  modal = (*inv_vdm) * nodal;
105  }
106 
107  // Operate with transpose of NodalToModal transformation
109  const Array<OneD, const NekDouble>& inarray,
110  Array<OneD, NekDouble>& outarray)
111  {
115  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
116 
117  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eCopy);
118  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
119  modal = Transpose(*inv_vdm) * nodal;
120  }
121 
123  const Array<OneD, const NekDouble>& inarray,
124  Array<OneD, NekDouble>& outarray)
125  {
126  StdMatrixKey Nkey(eNBasisTrans, DetShapeType(), *this,
129  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
130 
131  // Multiply out matrix
132  NekVector<NekDouble> modal(m_ncoeffs,inarray,eWrapper);
133  NekVector<NekDouble> nodal(m_ncoeffs,outarray,eWrapper);
134  nodal = (*vdm)*modal;
135  }
136 
141  {
142  LibUtilities::PointsManager()[m_nodalPointsKey]->GetPoints(x,y,z);
143  }
144 
146  {
147  int i,j;
150  DNekMatSharedPtr Mat;
151 
154  GetNodalPoints(r,s,t);
155 
156  //Store the values of m_phys in a temporary array
157  int nqtot = GetTotPoints();
158  Array<OneD,NekDouble> tmp_phys(nqtot);
159 
160  for(i = 0; i < m_ncoeffs; ++i)
161  {
162  // fill physical space with mode i
163  StdTetExp::v_FillMode(i,tmp_phys);
164 
165  // interpolate mode i to the Nodal points 'j' and
166  // store in outarray
167  for(j = 0; j < m_ncoeffs; ++j)
168  {
169  c[0] = r[j];
170  c[1] = s[j];
171  c[2] = t[j];
172  (*Mat)(j,i) = StdTetExp::v_PhysEvaluate(c,tmp_phys);
173  }
174  }
175 
176  return Mat;
177  }
178 
179 
180  //---------------------------------------
181  // Transforms
182  //---------------------------------------
183 
185  const Array<OneD, const NekDouble>& inarray,
186  Array<OneD, NekDouble>& outarray)
187  {
188  v_BwdTrans_SumFac(inarray,outarray);
189  }
190 
192  const Array<OneD, const NekDouble>& inarray,
193  Array<OneD, NekDouble>& outarray)
194  {
196  NodalToModal(inarray,tmp);
197  StdTetExp::v_BwdTrans_SumFac(tmp,outarray);
198  }
199 
201  const Array<OneD, const NekDouble>& inarray,
202  Array<OneD, NekDouble>& outarray)
203  {
204  v_IProductWRTBase(inarray,outarray);
205 
206  // get Mass matrix inverse
207  StdMatrixKey masskey(eInvMass, DetShapeType(), *this,
210  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
211 
212  // copy inarray in case inarray == outarray
213  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
215 
216  out = (*matsys)*in;
217  }
218 
219 
220  //---------------------------------------
221  // Inner product functions
222  //---------------------------------------
223 
225  const Array<OneD, const NekDouble>& inarray,
226  Array<OneD, NekDouble>& outarray)
227  {
228  v_IProductWRTBase_SumFac(inarray,outarray);
229  }
230 
232  const Array<OneD, const NekDouble>& inarray,
233  Array<OneD, NekDouble>& outarray,
234  bool multiplybyweights)
235  {
236  StdTetExp::v_IProductWRTBase_SumFac(inarray,outarray,multiplybyweights);
237  NodalToModalTranspose(outarray,outarray);
238  }
239 
241  const int dir,
242  const Array<OneD, const NekDouble>& inarray,
243  Array<OneD, NekDouble>& outarray)
244  {
245  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
246  }
247 
249  const int dir,
250  const Array<OneD, const NekDouble>& inarray,
251  Array<OneD, NekDouble>& outarray)
252  {
253  StdTetExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
254  NodalToModalTranspose(outarray,outarray);
255  }
256 
257  //---------------------------------------
258  // Evaluation functions
259  //---------------------------------------
260 
262  const int mode,
263  Array<OneD, NekDouble> &outarray)
264  {
265  ASSERTL2(mode >= m_ncoeffs,
266  "calling argument mode is larger than total expansion order");
267 
268  Vmath::Zero(m_ncoeffs, outarray, 1);
269  outarray[mode] = 1.0;
270  v_BwdTrans(outarray,outarray);
271  }
272 
273 
274  //---------------------------------------
275  // Mapping functions
276  //---------------------------------------
277 
278  /*
279  void StdNodalTriExp::v_GetFaceToElementMap(
280  const int fid,
281  const FaceOrientation faceOrient,
282  Array<OneD, unsigned int> &maparray,
283  Array<OneD, int> &signarray,
284  int nummodesA,
285  int nummodesB)
286  {
287  int P, Q, i, j, k, idx = 0, nFaceCoeffs = 0;
288 
289  ASSERTL0(fid >= 0 && fid <= 3,
290  "Local face ID must be between 0 and 3");
291 
292  if (nummodesA == -1)
293  {
294  switch(fid)
295  {
296  case 0:
297  nummodesA = m_base[0]->GetNumModes();
298  nummodesB = m_base[1]->GetNumModes();
299  break;
300  case 1:
301  nummodesA = m_base[0]->GetNumModes();
302  nummodesB = m_base[2]->GetNumModes();
303  break;
304  case 2:
305  case 3:
306  nummodesA = m_base[1]->GetNumModes();
307  nummodesB = m_base[2]->GetNumModes();
308  break;
309  }
310  }
311 
312  P = nummodesA;
313  Q = nummodesB;
314  nFaceCoeffs = Q + ((P-1)*(1 + 2*(Q-1) - (P-1)))/2;
315 
316  if (maparray.num_elements() != nFaceCoeffs)
317  {
318  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
319  }
320 
321  if (signarray.num_elements() != nFaceCoeffs)
322  {
323  signarray = Array<OneD, int>(nFaceCoeffs,1);
324  }
325  else
326  {
327  fill(signarray.get(), signarray.get()+nFaceCoeffs, 1);
328  }
329 
330  switch(fid)
331  {
332  case 0:
333  // Add vertices.
334  maparray[idx++] = 0;
335  maparray[idx++] = 1;
336  maparray[idx++] = 2;
337 
338  // Add edges.
339  for (i = 2; i < P; ++i)
340  {
341  maparray[idx++] = ;
342  }
343  }
344  }
345  */
346 
347  int StdNodalTetExp::v_GetVertexMap(const int localVertexId,
348  bool useCoeffPacking)
349  {
350  ASSERTL0(localVertexId >= 0 && localVertexId <= 3,
351  "Local Vertex ID must be between 0 and 3");
352  return localVertexId;
353  }
354 
356  Array<OneD, unsigned int>& outarray)
357  {
358  unsigned int i;
359  const unsigned int nBndryCoeff = NumBndryCoeffs();
360 
361  if (outarray.num_elements() != nBndryCoeff)
362  {
363  outarray = Array<OneD, unsigned int>(nBndryCoeff);
364  }
365 
366  for (i = 0; i < nBndryCoeff; i++)
367  {
368  outarray[i] = i;
369  }
370  }
371 
373  Array<OneD, unsigned int>& outarray)
374  {
375  unsigned int i;
376  const unsigned int nBndryCoeff = NumBndryCoeffs();
377 
378  if (outarray.num_elements() != m_ncoeffs-nBndryCoeff)
379  {
380  outarray = Array<OneD, unsigned int>(
381  m_ncoeffs-nBndryCoeff);
382  }
383 
384  for (i = nBndryCoeff; i < m_ncoeffs; i++)
385  {
386  outarray[i-nBndryCoeff] = i;
387  }
388  }
389 
390 
391  //---------------------------------------
392  // Wrapper functions
393  //---------------------------------------
394 
396  {
397  DNekMatSharedPtr Mat;
398 
399  switch(mkey.GetMatrixType())
400  {
401  case eNBasisTrans:
402  Mat = GenNBasisTransMatrix();
403  break;
404  default:
406  break;
407  }
408 
409  return Mat;
410  }
411 
413  const StdMatrixKey &mkey)
414  {
415  return StdNodalTetExp::v_GenMatrix(mkey);
416  }
417  } // end of namespace
418 } // end of namespace
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
void GetNodalPoints(Array< OneD, const NekDouble > &x, Array< OneD, const NekDouble > &y, Array< OneD, const NekDouble > &z)
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::PointsKey m_nodalPointsKey
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTetExp.cpp:545
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:70
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
virtual int v_GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:898
The base class for all shapes.
Definition: StdExpansion.h:69
void NodalToModal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:323
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:713
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
PointsManagerT & PointsManager(void)
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_BwdTrans(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 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.
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool mult=true)
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Describes the specification for a Basis.
Definition: Basis.h:50
PointsType GetPointsType() const
Definition: Points.h:111
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:226
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:249