Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 using namespace std;
40 
41 namespace Nektar
42 {
43  namespace StdRegions
44  {
45  StdNodalTriExp::StdNodalTriExp():
46  StdTriExp(),
47  m_nodalPointsKey()
48  {
49  }
50 
52  const LibUtilities::BasisKey &Ba,
53  const LibUtilities::BasisKey &Bb,
55  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(
56  Ba.GetNumModes(),
57  Bb.GetNumModes()),
58  2,Ba,Bb),
59  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
60  Ba.GetNumModes(),
61  Bb.GetNumModes()),
62  Ba,Bb),
63  StdTriExp (Ba,Bb),
64  m_nodalPointsKey(Ba.GetNumModes(),Ntype)
65  {
66  ASSERTL0(m_base[0]->GetNumModes() == m_base[1]->GetNumModes(),
67  "Nodal basis initiated with different orders in the a "
68  "and b directions");
69  }
70 
72  StdExpansion(T),
73  StdExpansion2D(T),
74  StdTriExp(T),
75  m_nodalPointsKey(T.m_nodalPointsKey)
76  {
77  }
78 
80  {
81  }
82 
84  {
85  return true;
86  }
87 
88  //-------------------------------
89  // Nodal basis specific routines
90  //-------------------------------
91 
93  const Array<OneD, const NekDouble>& inarray,
94  Array<OneD, NekDouble>& outarray)
95  {
99  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
100 
101  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eWrapper);
102  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
103  modal = (*inv_vdm) * nodal;
104  }
105 
106  // Operate with transpose of NodalToModal transformation
108  const Array<OneD, const NekDouble>& inarray,
109  Array<OneD, NekDouble>& outarray)
110  {
114  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
115 
116  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eCopy);
117  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
118  modal = Transpose(*inv_vdm) * nodal;
119  }
120 
122  const Array<OneD, const NekDouble>& inarray,
123  Array<OneD, NekDouble>& outarray)
124  {
125  StdMatrixKey Nkey(eNBasisTrans, DetShapeType(), *this,
128  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
129 
130  // Multiply out matrix
131  NekVector<NekDouble> modal(m_ncoeffs,inarray,eWrapper);
132  NekVector<NekDouble> nodal(m_ncoeffs,outarray,eWrapper);
133  nodal = (*vdm)*modal;
134  }
135 
139  {
141  }
142 
144  {
145  int i,j;
148  DNekMatSharedPtr Mat;
149 
152  GetNodalPoints(r,s);
153 
154  //Store the values of m_phys in a temporary array
155  int nqtot = GetTotPoints();
156  Array<OneD,NekDouble> phys(nqtot);
157 
158  for(i = 0; i < m_ncoeffs; ++i)
159  {
160  // fill physical space with mode i
161  StdTriExp::v_FillMode(i,phys);
162 
163  // interpolate mode i to the Nodal points 'j' and
164  // store in outarray
165  for(j = 0; j < m_ncoeffs; ++j)
166  {
167  c[0] = r[j];
168  c[1] = s[j];
169  (*Mat)(j,i) = StdTriExp::v_PhysEvaluate(c,phys);
170  }
171  }
172  return Mat;
173  }
174 
175 
176  //---------------------------------------
177  // Transforms
178  //---------------------------------------
179 
181  const Array<OneD, const NekDouble>& inarray,
182  Array<OneD, NekDouble>& outarray)
183  {
184  v_BwdTrans_SumFac(inarray,outarray);
185  }
186 
188  const Array<OneD, const NekDouble>& inarray,
189  Array<OneD, NekDouble>& outarray)
190  {
192  NodalToModal(inarray,tmp);
193  StdTriExp::v_BwdTrans_SumFac(tmp,outarray);
194  }
195 
197  const Array<OneD, const NekDouble>& inarray,
198  Array<OneD, NekDouble>& outarray)
199  {
200  v_IProductWRTBase(inarray,outarray);
201 
202  // get Mass matrix inverse
203  StdMatrixKey masskey(eInvMass, DetShapeType(), *this,
206  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
207 
208  // copy inarray in case inarray == outarray
209  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
211 
212  out = (*matsys)*in;
213  }
214 
215 
216  //---------------------------------------
217  // Inner product functions
218  //---------------------------------------
219 
221  const Array<OneD, const NekDouble>& inarray,
222  Array<OneD, NekDouble>& outarray)
223  {
224  v_IProductWRTBase_SumFac(inarray,outarray);
225  }
226 
228  const Array<OneD, const NekDouble>& inarray,
229  Array<OneD, NekDouble>& outarray,
230  bool multiplybyweights)
231  {
232  StdTriExp::v_IProductWRTBase_SumFac(inarray,outarray,multiplybyweights);
233  NodalToModalTranspose(outarray,outarray);
234  }
235 
237  const int dir,
238  const Array<OneD, const NekDouble>& inarray,
239  Array<OneD, NekDouble>& outarray)
240  {
241  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
242  }
243 
245  const int dir,
246  const Array<OneD, const NekDouble>& inarray,
247  Array<OneD, NekDouble>& outarray)
248  {
249  StdTriExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
250  NodalToModalTranspose(outarray,outarray);
251  }
252 
253  //---------------------------------------
254  // Evaluation functions
255  //---------------------------------------
256 
258  const int mode,
259  Array<OneD, NekDouble> &outarray)
260  {
261  ASSERTL2(mode >= m_ncoeffs,
262  "calling argument mode is larger than total expansion order");
263 
264  Vmath::Zero(m_ncoeffs, outarray, 1);
265  outarray[mode] = 1.0;
266  v_BwdTrans(outarray,outarray);
267  }
268 
269 
270  //---------------------------
271  // Helper functions
272  //---------------------------
273 
275  {
276  return 3 + (GetBasisNumModes(0)-2) + 2*(GetBasisNumModes(1)-2);
277  }
278 
279  //--------------------------
280  // Mappings
281  //--------------------------
282 
284  const int eid,
285  const Orientation edgeOrient,
286  Array<OneD, unsigned int> &maparray,
287  Array<OneD, int> &signarray,
288  int P)
289  {
290  ASSERTL0(eid >= 0 && eid <= 2,
291  "Local Edge ID must be between 0 and 2");
292 
293  ASSERTL0(P == -1, "Nodal triangle not set up to deal with variable"
294  "polynomial order.");
295  const int nEdgeCoeffs = GetEdgeNcoeffs(eid);
296 
297  if (maparray.num_elements() != nEdgeCoeffs)
298  {
299  maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
300  }
301 
302  if (signarray.num_elements() != nEdgeCoeffs)
303  {
304  signarray = Array<OneD, int>(nEdgeCoeffs,1);
305  }
306  else
307  {
308  fill(signarray.get(), signarray.get()+nEdgeCoeffs, 1);
309  }
310 
311  maparray[0] = eid;
312  maparray[1] = eid == 2 ? 0 : eid+1;
313  for (int i = 2; i < nEdgeCoeffs; i++)
314  {
315  maparray[i] = eid*(nEdgeCoeffs-2)+1+i;
316  }
317 
318  if (edgeOrient == eBackwards)
319  {
320  reverse(maparray.get(), maparray.get()+nEdgeCoeffs);
321  }
322  }
323 
324  int StdNodalTriExp::v_GetVertexMap(const int localVertexId,
325  bool useCoeffPacking)
326  {
327  ASSERTL0(localVertexId >= 0 && localVertexId <= 2,
328  "Local Vertex ID must be between 0 and 2");
329  return localVertexId;
330  }
331 
333  const int eid,
334  const Orientation edgeOrient,
335  Array<OneD, unsigned int> &maparray,
336  Array<OneD, int> &signarray)
337  {
338  ASSERTL0(eid >= 0 && eid <= 2,
339  "Local Edge ID must be between 0 and 2");
340 
341  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
342 
343  if (maparray.num_elements() != nEdgeIntCoeffs)
344  {
345  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
346  }
347 
348  if (signarray.num_elements() != nEdgeIntCoeffs)
349  {
350  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
351  }
352  else
353  {
354  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
355  }
356 
357  for (int i = 0; i < nEdgeIntCoeffs; i++)
358  {
359  maparray[i] = eid*nEdgeIntCoeffs+3+i;
360  }
361 
362  if (edgeOrient == eBackwards)
363  {
364  reverse(maparray.get(), maparray.get()+nEdgeIntCoeffs);
365  }
366  }
367 
369  Array<OneD, unsigned int>& outarray)
370  {
371  unsigned int i;
372  if (outarray.num_elements() != GetNcoeffs()-NumBndryCoeffs())
373  {
374  outarray = Array<OneD, unsigned int>(
376  }
377 
378  for (i = NumBndryCoeffs(); i < GetNcoeffs(); i++)
379  {
380  outarray[i-NumBndryCoeffs()] = i;
381  }
382  }
383 
385  Array<OneD, unsigned int>& outarray)
386  {
387  unsigned int i;
388  if (outarray.num_elements()!=NumBndryCoeffs())
389  {
391  }
392 
393  for (i = 0; i < NumBndryCoeffs(); i++)
394  {
395  outarray[i] = i;
396  }
397  }
398 
399 
400  //---------------------------------------
401  // Wrapper functions
402  //---------------------------------------
403 
405  {
406  DNekMatSharedPtr Mat;
407 
408  switch(mkey.GetMatrixType())
409  {
410  case eNBasisTrans:
411  Mat = GenNBasisTransMatrix();
412  break;
413  default:
415  break;
416  }
417 
418  return Mat;
419  }
420 
422  const StdMatrixKey &mkey)
423  {
424  return StdNodalTriExp::v_GenMatrix(mkey);
425  }
426 
427 
428  //---------------------------------------
429  // Operator evaluation functions
430  //---------------------------------------
431 
433  const Array<OneD, const NekDouble> &inarray,
434  Array<OneD, NekDouble> &outarray,
435  const StdMatrixKey &mkey)
436  {
437  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
438  }
439 
441  const Array<OneD, const NekDouble> &inarray,
442  Array<OneD, NekDouble> &outarray,
443  const StdMatrixKey &mkey)
444  {
446  inarray,outarray,mkey);
447  }
448 
450  const int k1,
451  const int k2,
452  const Array<OneD, const NekDouble> &inarray,
453  Array<OneD, NekDouble> &outarray,
454  const StdMatrixKey &mkey)
455 
456  {
458  k1,k2,inarray,outarray,mkey);
459  }
460 
462  const int i,
463  const Array<OneD, const NekDouble> &inarray,
464  Array<OneD, NekDouble> &outarray,
465  const StdMatrixKey &mkey)
466  {
467  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
468  }
469 
471  const Array<OneD, const NekDouble> &inarray,
472  Array<OneD, NekDouble> &outarray,
473  const StdMatrixKey &mkey)
474  {
476  inarray,outarray,mkey);
477  }
478 
479 
480  //---------------------------------------
481  // Private helper functions
482  //---------------------------------------
483 
484  } // end of namespace
485 } // end of namespace
486 
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:198
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)
STL namespace.
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:589
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:706
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:479
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:254
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:250
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:686
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:373
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:228
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:253