Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdNodalTriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdNodalTriExp.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 triangle routines built upon StdExpansion2D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 
39 namespace Nektar
40 {
41  namespace StdRegions
42  {
44  StdTriExp(),
45  m_nodalPointsKey()
46  {
47  }
48 
50  const LibUtilities::BasisKey &Ba,
51  const LibUtilities::BasisKey &Bb,
53  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(
54  Ba.GetNumModes(),
55  Bb.GetNumModes()),
56  2,Ba,Bb),
57  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
58  Ba.GetNumModes(),
59  Bb.GetNumModes()),
60  Ba,Bb),
61  StdTriExp (Ba,Bb),
62  m_nodalPointsKey(Ba.GetNumModes(),Ntype)
63  {
64  ASSERTL0(m_base[0]->GetNumModes() == m_base[1]->GetNumModes(),
65  "Nodal basis initiated with different orders in the a "
66  "and b directions");
67  }
68 
70  StdExpansion(T),
71  StdExpansion2D(T),
72  StdTriExp(T),
73  m_nodalPointsKey(T.m_nodalPointsKey)
74  {
75  }
76 
78  {
79  }
80 
82  {
83  return true;
84  }
85 
86  //-------------------------------
87  // Nodal basis specific routines
88  //-------------------------------
89 
91  const Array<OneD, const NekDouble>& inarray,
92  Array<OneD, NekDouble>& outarray)
93  {
97  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
98 
100  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
101  modal = (*inv_vdm) * nodal;
102  }
103 
104  // Operate with transpose of NodalToModal transformation
106  const Array<OneD, const NekDouble>& inarray,
107  Array<OneD, NekDouble>& outarray)
108  {
112  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
113 
114  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eCopy);
115  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
116  modal = Transpose(*inv_vdm) * nodal;
117  }
118 
120  const Array<OneD, const NekDouble>& inarray,
121  Array<OneD, NekDouble>& outarray)
122  {
123  StdMatrixKey Nkey(eNBasisTrans, DetShapeType(), *this,
126  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
127 
128  // Multiply out matrix
129  NekVector<NekDouble> modal(m_ncoeffs,inarray,eWrapper);
130  NekVector<NekDouble> nodal(m_ncoeffs,outarray,eWrapper);
131  nodal = (*vdm)*modal;
132  }
133 
137  {
139  }
140 
142  {
143  int i,j;
146  DNekMatSharedPtr Mat;
147 
150  GetNodalPoints(r,s);
151 
152  //Store the values of m_phys in a temporary array
153  int nqtot = GetTotPoints();
154  Array<OneD,NekDouble> phys(nqtot);
155 
156  for(i = 0; i < m_ncoeffs; ++i)
157  {
158  // fill physical space with mode i
159  StdTriExp::v_FillMode(i,phys);
160 
161  // interpolate mode i to the Nodal points 'j' and
162  // store in outarray
163  for(j = 0; j < m_ncoeffs; ++j)
164  {
165  c[0] = r[j];
166  c[1] = s[j];
167  (*Mat)(j,i) = StdTriExp::v_PhysEvaluate(c,phys);
168  }
169  }
170  return Mat;
171  }
172 
173 
174  //---------------------------------------
175  // Transforms
176  //---------------------------------------
177 
179  const Array<OneD, const NekDouble>& inarray,
180  Array<OneD, NekDouble>& outarray)
181  {
182  v_BwdTrans_SumFac(inarray,outarray);
183  }
184 
186  const Array<OneD, const NekDouble>& inarray,
187  Array<OneD, NekDouble>& outarray)
188  {
190  NodalToModal(inarray,tmp);
191  StdTriExp::v_BwdTrans_SumFac(tmp,outarray);
192  }
193 
195  const Array<OneD, const NekDouble>& inarray,
196  Array<OneD, NekDouble>& outarray)
197  {
198  v_IProductWRTBase(inarray,outarray);
199 
200  // get Mass matrix inverse
201  StdMatrixKey masskey(eInvMass, DetShapeType(), *this,
204  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
205 
206  // copy inarray in case inarray == outarray
207  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
209 
210  out = (*matsys)*in;
211  }
212 
213 
214  //---------------------------------------
215  // Inner product functions
216  //---------------------------------------
217 
219  const Array<OneD, const NekDouble>& inarray,
220  Array<OneD, NekDouble>& outarray)
221  {
222  v_IProductWRTBase_SumFac(inarray,outarray);
223  }
224 
226  const Array<OneD, const NekDouble>& inarray,
227  Array<OneD, NekDouble>& outarray,
228  bool multiplybyweights)
229  {
230  StdTriExp::v_IProductWRTBase_SumFac(inarray,outarray,multiplybyweights);
231  NodalToModalTranspose(outarray,outarray);
232  }
233 
235  const int dir,
236  const Array<OneD, const NekDouble>& inarray,
237  Array<OneD, NekDouble>& outarray)
238  {
239  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
240  }
241 
243  const int dir,
244  const Array<OneD, const NekDouble>& inarray,
245  Array<OneD, NekDouble>& outarray)
246  {
247  StdTriExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
248  NodalToModalTranspose(outarray,outarray);
249  }
250 
251  //---------------------------------------
252  // Evaluation functions
253  //---------------------------------------
254 
256  const int mode,
257  Array<OneD, NekDouble> &outarray)
258  {
259  ASSERTL2(mode >= m_ncoeffs,
260  "calling argument mode is larger than total expansion order");
261 
262  Vmath::Zero(m_ncoeffs, outarray, 1);
263  outarray[mode] = 1.0;
264  v_BwdTrans(outarray,outarray);
265  }
266 
267 
268  //---------------------------
269  // Helper functions
270  //---------------------------
271 
273  {
274  return 3 + (GetBasisNumModes(0)-2) + 2*(GetBasisNumModes(1)-2);
275  }
276 
277  //--------------------------
278  // Mappings
279  //--------------------------
280 
282  const int eid,
283  const Orientation edgeOrient,
284  Array<OneD, unsigned int> &maparray,
285  Array<OneD, int> &signarray,
286  int P)
287  {
288  ASSERTL0(eid >= 0 && eid <= 2,
289  "Local Edge ID must be between 0 and 2");
290 
291  ASSERTL0(P == -1, "Nodal triangle not set up to deal with variable"
292  "polynomial order.");
293  const int nEdgeCoeffs = GetEdgeNcoeffs(eid);
294 
295  if (maparray.num_elements() != nEdgeCoeffs)
296  {
297  maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
298  }
299 
300  if (signarray.num_elements() != nEdgeCoeffs)
301  {
302  signarray = Array<OneD, int>(nEdgeCoeffs,1);
303  }
304  else
305  {
306  fill(signarray.get(), signarray.get()+nEdgeCoeffs, 1);
307  }
308 
309  maparray[0] = eid;
310  maparray[1] = eid == 2 ? 0 : eid+1;
311  for (int i = 2; i < nEdgeCoeffs; i++)
312  {
313  maparray[i] = eid*(nEdgeCoeffs-2)+1+i;
314  }
315 
316  if (edgeOrient == eBackwards)
317  {
318  reverse(maparray.get(), maparray.get()+nEdgeCoeffs);
319  }
320  }
321 
322  int StdNodalTriExp::v_GetVertexMap(const int localVertexId,
323  bool useCoeffPacking)
324  {
325  ASSERTL0(localVertexId >= 0 && localVertexId <= 2,
326  "Local Vertex ID must be between 0 and 2");
327  return localVertexId;
328  }
329 
331  const int eid,
332  const Orientation edgeOrient,
333  Array<OneD, unsigned int> &maparray,
334  Array<OneD, int> &signarray)
335  {
336  ASSERTL0(eid >= 0 && eid <= 2,
337  "Local Edge ID must be between 0 and 2");
338 
339  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
340 
341  if (maparray.num_elements() != nEdgeIntCoeffs)
342  {
343  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
344  }
345 
346  if (signarray.num_elements() != nEdgeIntCoeffs)
347  {
348  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
349  }
350  else
351  {
352  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
353  }
354 
355  for (int i = 0; i < nEdgeIntCoeffs; i++)
356  {
357  maparray[i] = eid*nEdgeIntCoeffs+3+i;
358  }
359 
360  if (edgeOrient == eBackwards)
361  {
362  reverse(maparray.get(), maparray.get()+nEdgeIntCoeffs);
363  }
364  }
365 
367  Array<OneD, unsigned int>& outarray)
368  {
369  unsigned int i;
370  if (outarray.num_elements() != GetNcoeffs()-NumBndryCoeffs())
371  {
372  outarray = Array<OneD, unsigned int>(
374  }
375 
376  for (i = NumBndryCoeffs(); i < GetNcoeffs(); i++)
377  {
378  outarray[i-NumBndryCoeffs()] = i;
379  }
380  }
381 
383  Array<OneD, unsigned int>& outarray)
384  {
385  unsigned int i;
386  if (outarray.num_elements()!=NumBndryCoeffs())
387  {
389  }
390 
391  for (i = 0; i < NumBndryCoeffs(); i++)
392  {
393  outarray[i] = i;
394  }
395  }
396 
397 
398  //---------------------------------------
399  // Wrapper functions
400  //---------------------------------------
401 
403  {
404  DNekMatSharedPtr Mat;
405 
406  switch(mkey.GetMatrixType())
407  {
408  case eNBasisTrans:
409  Mat = GenNBasisTransMatrix();
410  break;
411  default:
413  break;
414  }
415 
416  return Mat;
417  }
418 
420  const StdMatrixKey &mkey)
421  {
422  return StdNodalTriExp::v_GenMatrix(mkey);
423  }
424 
425 
426  //---------------------------------------
427  // Operator evaluation functions
428  //---------------------------------------
429 
431  const Array<OneD, const NekDouble> &inarray,
432  Array<OneD, NekDouble> &outarray,
433  const StdMatrixKey &mkey)
434  {
435  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
436  }
437 
439  const Array<OneD, const NekDouble> &inarray,
440  Array<OneD, NekDouble> &outarray,
441  const StdMatrixKey &mkey)
442  {
444  inarray,outarray,mkey);
445  }
446 
448  const int k1,
449  const int k2,
450  const Array<OneD, const NekDouble> &inarray,
451  Array<OneD, NekDouble> &outarray,
452  const StdMatrixKey &mkey)
453 
454  {
456  k1,k2,inarray,outarray,mkey);
457  }
458 
460  const int i,
461  const Array<OneD, const NekDouble> &inarray,
462  Array<OneD, NekDouble> &outarray,
463  const StdMatrixKey &mkey)
464  {
465  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
466  }
467 
469  const Array<OneD, const NekDouble> &inarray,
470  Array<OneD, NekDouble> &outarray,
471  const StdMatrixKey &mkey)
472  {
474  inarray,outarray,mkey);
475  }
476 
477 
478  //---------------------------------------
479  // Private helper functions
480  //---------------------------------------
481 
482  } // end of namespace
483 } // end of namespace
484 
void HelmholtzMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void LaplacianMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &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[p]*base1[pq] and put into ou...
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:178
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:593
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion.h:287
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:700
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward tranform for triangular elements.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
LibUtilities::PointsKey m_nodalPointsKey
The base class for all shapes.
Definition: StdExpansion.h:69
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
virtual void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Transform a given function from physical quadrature space to coefficient space.
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void GetNodalPoints(Array< OneD, const NekDouble > &x, Array< OneD, const NekDouble > &y)
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.
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
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_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTriExp.cpp:483
virtual void v_IProductWRTDerivBase_SumFac(const int dir, 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: StdTriExp.cpp:258
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:690
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void NodalToModal(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
Array< OneD, LibUtilities::BasisSharedPtr > m_base
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Describes the specification for a Basis.
Definition: Basis.h:50
PointsType GetPointsType() const
Definition: Points.h:111
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:227
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:252