Nektar++
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 // 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 triangle routines built upon StdExpansion2D
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 namespace StdRegions
45 {
46 StdNodalTriExp::StdNodalTriExp() : StdTriExp(), m_nodalPointsKey()
47 {
48 }
49 
51  const LibUtilities::BasisKey &Bb,
53  : StdExpansion(LibUtilities::StdTriData::getNumberOfCoefficients(
54  Ba.GetNumModes(), Bb.GetNumModes()),
55  2, Ba, Bb),
56  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
57  Ba.GetNumModes(), Bb.GetNumModes()),
58  Ba, Bb),
59  StdTriExp(Ba, Bb), m_nodalPointsKey(Ba.GetNumModes(), Ntype)
60 {
61  ASSERTL0(m_base[0]->GetNumModes() == m_base[1]->GetNumModes(),
62  "Nodal basis initiated with different orders in the a "
63  "and b directions");
64 }
65 
68  m_nodalPointsKey(T.m_nodalPointsKey)
69 {
70 }
71 
73 {
74 }
75 
77 {
78  return true;
79 }
80 
81 //-------------------------------
82 // Nodal basis specific routines
83 //-------------------------------
84 
86  Array<OneD, NekDouble> &outarray)
87 {
91  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
92 
93  NekVector<NekDouble> nodal(m_ncoeffs, inarray, eWrapper);
94  NekVector<NekDouble> modal(m_ncoeffs, outarray, eWrapper);
95  modal = (*inv_vdm) * nodal;
96 }
97 
98 // Operate with transpose of NodalToModal transformation
100  const Array<OneD, const NekDouble> &inarray,
101  Array<OneD, NekDouble> &outarray)
102 {
106  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
107 
108  NekVector<NekDouble> nodal(m_ncoeffs, inarray, eCopy);
109  NekVector<NekDouble> modal(m_ncoeffs, outarray, eWrapper);
110  modal = Transpose(*inv_vdm) * nodal;
111 }
112 
114  Array<OneD, NekDouble> &outarray)
115 {
118  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
119 
120  // Multiply out matrix
121  NekVector<NekDouble> modal(m_ncoeffs, inarray, eWrapper);
122  NekVector<NekDouble> nodal(m_ncoeffs, outarray, eWrapper);
123  nodal = (*vdm) * modal;
124 }
125 
128 {
129  LibUtilities::PointsManager()[m_nodalPointsKey]->GetPoints(x, y);
130 }
131 
133 {
134  int i, j;
137  DNekMatSharedPtr Mat;
138 
140  GetNodalPoints(r, s);
141 
142  // Store the values of m_phys in a temporary array
143  int nqtot = GetTotPoints();
144  Array<OneD, NekDouble> phys(nqtot);
145 
146  for (i = 0; i < m_ncoeffs; ++i)
147  {
148  // fill physical space with mode i
149  StdTriExp::v_FillMode(i, phys);
150 
151  // interpolate mode i to the Nodal points 'j' and
152  // store in outarray
153  for (j = 0; j < m_ncoeffs; ++j)
154  {
155  c[0] = r[j];
156  c[1] = s[j];
157  (*Mat)(j, i) = StdTriExp::v_PhysEvaluate(c, phys);
158  }
159  }
160  return Mat;
161 }
162 
163 //---------------------------------------
164 // Transforms
165 //---------------------------------------
166 
168  Array<OneD, NekDouble> &outarray)
169 {
170  v_BwdTrans_SumFac(inarray, outarray);
171 }
172 
174  const Array<OneD, const NekDouble> &inarray,
175  Array<OneD, NekDouble> &outarray)
176 {
178  NodalToModal(inarray, tmp);
179  StdTriExp::v_BwdTrans_SumFac(tmp, outarray);
180 }
181 
183  Array<OneD, NekDouble> &outarray)
184 {
185  v_IProductWRTBase(inarray, outarray);
186 
187  // get Mass matrix inverse
190  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
191 
192  // copy inarray in case inarray == outarray
193  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
194  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
195 
196  out = (*matsys) * in;
197 }
198 
199 //---------------------------------------
200 // Inner product functions
201 //---------------------------------------
202 
204  const Array<OneD, const NekDouble> &inarray,
205  Array<OneD, NekDouble> &outarray)
206 {
207  v_IProductWRTBase_SumFac(inarray, outarray);
208 }
209 
211  const Array<OneD, const NekDouble> &inarray,
212  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
213 {
214  StdTriExp::v_IProductWRTBase_SumFac(inarray, outarray, multiplybyweights);
215  NodalToModalTranspose(outarray, outarray);
216 }
217 
219  const int dir, const Array<OneD, const NekDouble> &inarray,
220  Array<OneD, NekDouble> &outarray)
221 {
222  v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
223 }
224 
226  const int dir, const Array<OneD, const NekDouble> &inarray,
227  Array<OneD, NekDouble> &outarray)
228 {
229  StdTriExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
230  NodalToModalTranspose(outarray, outarray);
231 }
232 
233 //---------------------------------------
234 // Evaluation functions
235 //---------------------------------------
236 
237 void StdNodalTriExp::v_FillMode(const int mode,
238  Array<OneD, NekDouble> &outarray)
239 {
240  ASSERTL2(mode >= m_ncoeffs,
241  "calling argument mode is larger than total expansion order");
242 
243  Vmath::Zero(m_ncoeffs, outarray, 1);
244  outarray[mode] = 1.0;
245  v_BwdTrans(outarray, outarray);
246 }
247 
248 //---------------------------
249 // Helper functions
250 //---------------------------
251 
253 {
254  return 3 + (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
255 }
256 
257 //--------------------------
258 // Mappings
259 //--------------------------
260 
261 int StdNodalTriExp::v_GetVertexMap(const int localVertexId,
262  bool useCoeffPacking)
263 {
264  boost::ignore_unused(useCoeffPacking);
265  ASSERTL0(localVertexId >= 0 && localVertexId <= 2,
266  "Local Vertex ID must be between 0 and 2");
267  return localVertexId;
268 }
269 
271  Array<OneD, unsigned int> &maparray,
272  Array<OneD, int> &signarray,
273  Orientation edgeOrient, int P,
274  int Q)
275 {
276  boost::ignore_unused(Q);
277 
278  ASSERTL0(eid >= 0 && eid <= 2, "Local Edge ID must be between 0 and 2");
279 
280  const int nEdgeCoeffs = GetTraceNcoeffs(eid);
281 
282  ASSERTL0(P == -1 || P == nEdgeCoeffs,
283  "Nodal triangle not set up to deal with variable "
284  "polynomial order.");
285 
286  if (maparray.size() != nEdgeCoeffs)
287  {
288  maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
289  }
290 
291  if (signarray.size() != nEdgeCoeffs)
292  {
293  signarray = Array<OneD, int>(nEdgeCoeffs, 1);
294  }
295  else
296  {
297  fill(signarray.get(), signarray.get() + nEdgeCoeffs, 1);
298  }
299 
300  Orientation orient = edgeOrient;
301  if (eid == 2)
302  {
303  orient = orient == eForwards ? eBackwards : eForwards;
304  }
305 
306  maparray[0] = eid;
307  maparray[nEdgeCoeffs - 1] = eid == 2 ? 0 : eid + 1;
308  for (int i = 2; i < nEdgeCoeffs; i++)
309  {
310  maparray[i - 1] = eid * (nEdgeCoeffs - 2) + 1 + i;
311  }
312 
313  if (orient == eBackwards)
314  {
315  reverse(maparray.get(), maparray.get() + nEdgeCoeffs);
316  }
317 }
318 
320  const int eid, Array<OneD, unsigned int> &maparray,
321  Array<OneD, int> &signarray, const Orientation edgeOrient)
322 {
323  ASSERTL0(eid >= 0 && eid <= 2, "Local Edge ID must be between 0 and 2");
324 
325  const int nEdgeIntCoeffs = GetTraceNcoeffs(eid) - 2;
326 
327  if (maparray.size() != nEdgeIntCoeffs)
328  {
329  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
330  }
331 
332  if (signarray.size() != nEdgeIntCoeffs)
333  {
334  signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
335  }
336  else
337  {
338  fill(signarray.get(), signarray.get() + nEdgeIntCoeffs, 1);
339  }
340 
341  Orientation orient = edgeOrient;
342  if (eid == 2)
343  {
344  orient = orient == eForwards ? eBackwards : eForwards;
345  }
346 
347  for (int i = 0; i < nEdgeIntCoeffs; i++)
348  {
349  maparray[i] = eid * nEdgeIntCoeffs + 3 + i;
350  }
351 
352  if (orient == eBackwards)
353  {
354  reverse(maparray.get(), maparray.get() + nEdgeIntCoeffs);
355  }
356 }
357 
359 {
360  unsigned int i;
361  if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
362  {
364  }
365 
366  for (i = NumBndryCoeffs(); i < GetNcoeffs(); i++)
367  {
368  outarray[i - NumBndryCoeffs()] = i;
369  }
370 }
371 
373 {
374  unsigned int i;
375  if (outarray.size() != NumBndryCoeffs())
376  {
378  }
379 
380  for (i = 0; i < NumBndryCoeffs(); i++)
381  {
382  outarray[i] = i;
383  }
384 }
385 
386 //---------------------------------------
387 // Wrapper functions
388 //---------------------------------------
389 
391 {
392  DNekMatSharedPtr Mat;
393 
394  switch (mkey.GetMatrixType())
395  {
396  case eNBasisTrans:
397  Mat = GenNBasisTransMatrix();
398  break;
399  default:
401  break;
402  }
403 
404  return Mat;
405 }
406 
408 {
409  return StdNodalTriExp::v_GenMatrix(mkey);
410 }
411 
412 //---------------------------------------
413 // Operator evaluation functions
414 //---------------------------------------
415 
417  Array<OneD, NekDouble> &outarray,
418  const StdMatrixKey &mkey)
419 {
420  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
421 }
422 
424  const Array<OneD, const NekDouble> &inarray,
425  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
426 {
428  mkey);
429 }
430 
432  const int k1, const int k2, const Array<OneD, const NekDouble> &inarray,
433  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
434 
435 {
436  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
437 }
438 
440  const int i, const Array<OneD, const NekDouble> &inarray,
441  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
442 {
443  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
444 }
445 
447  const Array<OneD, const NekDouble> &inarray,
448  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
449 {
451  mkey);
452 }
453 
454 //---------------------------------------
455 // Private helper functions
456 //---------------------------------------
457 
458 } // namespace StdRegions
459 } // 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
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 GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void LaplacianMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
void LaplacianMatrixOp_MatFree(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:375
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:269
void HelmholtzMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:176
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void NodalToModal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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)
virtual void v_LaplacianMatrixOp(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)
void GetNodalPoints(Array< OneD, const NekDouble > &x, Array< OneD, const NekDouble > &y)
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::PointsKey m_nodalPointsKey
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
virtual void v_GetTraceToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation edgeOrient=eForwards, int P=-1, int Q=-1)
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[p]*base1[pq] and put into ou...
virtual void v_GetTraceInteriorToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation edgeOrient=eForwards)
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward tranform for triangular elements.
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:246
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTriExp.cpp:462
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:673
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:568
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