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 // 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 tetrahedral 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::StdTetData::getNumberOfCoefficients(
49  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
50  3, Ba, Bb, Bc),
51  StdExpansion3D(LibUtilities::StdTetData::getNumberOfCoefficients(
52  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
53  Ba, Bb, Bc),
54  StdTetExp(Ba, Bb, Bc), m_nodalPointsKey(Ba.GetNumModes(), Ntype)
55 {
56  ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
57  "order in 'a' direction is higher than order "
58  "in 'b' direction");
59  ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
60  "order in 'a' direction is higher than order "
61  "in 'c' direction");
62  ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(),
63  "order in 'b' direction is higher than order "
64  "in 'c' direction");
65 }
66 
69  m_nodalPointsKey(T.m_nodalPointsKey)
70 {
71 }
72 
74 {
75 }
76 
78 {
79  return true;
80 }
81 
82 //-------------------------------
83 // Nodal basis specific routines
84 //-------------------------------
85 
87  Array<OneD, NekDouble> &outarray)
88 {
92  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
93 
94  NekVector<NekDouble> nodal(m_ncoeffs, inarray, eWrapper);
95  NekVector<NekDouble> modal(m_ncoeffs, outarray, eWrapper);
96  modal = (*inv_vdm) * nodal;
97 }
98 
99 // Operate with transpose of NodalToModal transformation
101  const Array<OneD, const NekDouble> &inarray,
102  Array<OneD, NekDouble> &outarray)
103 {
107  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
108 
109  NekVector<NekDouble> nodal(m_ncoeffs, inarray, eCopy);
110  NekVector<NekDouble> modal(m_ncoeffs, outarray, eWrapper);
111  modal = Transpose(*inv_vdm) * nodal;
112 }
113 
115  Array<OneD, NekDouble> &outarray)
116 {
119  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
120 
121  // Multiply out matrix
122  NekVector<NekDouble> modal(m_ncoeffs, inarray, eWrapper);
123  NekVector<NekDouble> nodal(m_ncoeffs, outarray, eWrapper);
124  nodal = (*vdm) * modal;
125 }
126 
130 {
131  LibUtilities::PointsManager()[m_nodalPointsKey]->GetPoints(x, y, z);
132 }
133 
135 {
136  int i, j;
139  DNekMatSharedPtr Mat;
140 
142  GetNodalPoints(r, s, t);
143 
144  // Store the values of m_phys in a temporary array
145  int nqtot = GetTotPoints();
146  Array<OneD, NekDouble> tmp_phys(nqtot);
147 
148  for (i = 0; i < m_ncoeffs; ++i)
149  {
150  // fill physical space with mode i
151  StdTetExp::v_FillMode(i, tmp_phys);
152 
153  // interpolate mode i to the Nodal points 'j' and
154  // store in outarray
155  for (j = 0; j < m_ncoeffs; ++j)
156  {
157  c[0] = r[j];
158  c[1] = s[j];
159  c[2] = t[j];
160  (*Mat)(j, i) = StdTetExp::v_PhysEvaluate(c, tmp_phys);
161  }
162  }
163 
164  return Mat;
165 }
166 
167 //---------------------------------------
168 // Transforms
169 //---------------------------------------
170 
172  Array<OneD, NekDouble> &outarray)
173 {
174  v_BwdTrans_SumFac(inarray, outarray);
175 }
176 
178  const Array<OneD, const NekDouble> &inarray,
179  Array<OneD, NekDouble> &outarray)
180 {
182  NodalToModal(inarray, tmp);
183  StdTetExp::v_BwdTrans_SumFac(tmp, outarray);
184 }
185 
187  Array<OneD, NekDouble> &outarray)
188 {
189  v_IProductWRTBase(inarray, outarray);
190 
191  // get Mass matrix inverse
194  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
195 
196  // copy inarray in case inarray == outarray
197  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
198  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
199 
200  out = (*matsys) * in;
201 }
202 
203 //---------------------------------------
204 // Inner product functions
205 //---------------------------------------
206 
208  const Array<OneD, const NekDouble> &inarray,
209  Array<OneD, NekDouble> &outarray)
210 {
211  v_IProductWRTBase_SumFac(inarray, outarray);
212 }
213 
215  const Array<OneD, const NekDouble> &inarray,
216  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
217 {
218  StdTetExp::v_IProductWRTBase_SumFac(inarray, outarray, multiplybyweights);
219  NodalToModalTranspose(outarray, outarray);
220 }
221 
223  const int dir, const Array<OneD, const NekDouble> &inarray,
224  Array<OneD, NekDouble> &outarray)
225 {
226  v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
227 }
228 
230  const int dir, const Array<OneD, const NekDouble> &inarray,
231  Array<OneD, NekDouble> &outarray)
232 {
233  StdTetExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
234  NodalToModalTranspose(outarray, outarray);
235 }
236 
237 //---------------------------------------
238 // Evaluation functions
239 //---------------------------------------
240 
241 void StdNodalTetExp::v_FillMode(const int mode,
242  Array<OneD, NekDouble> &outarray)
243 {
244  ASSERTL2(mode >= m_ncoeffs,
245  "calling argument mode is larger than total expansion order");
246 
247  Vmath::Zero(m_ncoeffs, outarray, 1);
248  outarray[mode] = 1.0;
249  v_BwdTrans(outarray, outarray);
250 }
251 
252 //---------------------------------------
253 // Mapping functions
254 //---------------------------------------
255 
256 /*
257 void StdNodalTriExp::v_GetFaceToElementMap(
258  const int fid,
259  const FaceOrientation faceOrient,
260  Array<OneD, unsigned int> &maparray,
261  Array<OneD, int> &signarray,
262  int nummodesA,
263  int nummodesB)
264 {
265  int P, Q, i, j, k, idx = 0, nFaceCoeffs = 0;
266 
267  ASSERTL0(fid >= 0 && fid <= 3,
268  "Local face ID must be between 0 and 3");
269 
270  if (nummodesA == -1)
271  {
272  switch(fid)
273  {
274  case 0:
275  nummodesA = m_base[0]->GetNumModes();
276  nummodesB = m_base[1]->GetNumModes();
277  break;
278  case 1:
279  nummodesA = m_base[0]->GetNumModes();
280  nummodesB = m_base[2]->GetNumModes();
281  break;
282  case 2:
283  case 3:
284  nummodesA = m_base[1]->GetNumModes();
285  nummodesB = m_base[2]->GetNumModes();
286  break;
287  }
288  }
289 
290  P = nummodesA;
291  Q = nummodesB;
292  nFaceCoeffs = Q + ((P-1)*(1 + 2*(Q-1) - (P-1)))/2;
293 
294  if (maparray.size() != nFaceCoeffs)
295  {
296  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
297  }
298 
299  if (signarray.size() != nFaceCoeffs)
300  {
301  signarray = Array<OneD, int>(nFaceCoeffs,1);
302  }
303  else
304  {
305  fill(signarray.get(), signarray.get()+nFaceCoeffs, 1);
306  }
307 
308  switch(fid)
309  {
310  case 0:
311  // Add vertices.
312  maparray[idx++] = 0;
313  maparray[idx++] = 1;
314  maparray[idx++] = 2;
315 
316  // Add edges.
317  for (i = 2; i < P; ++i)
318  {
319  maparray[idx++] = ;
320  }
321  }
322 }
323 */
324 
325 int StdNodalTetExp::v_GetVertexMap(const int localVertexId,
326  bool useCoeffPacking)
327 {
328  boost::ignore_unused(useCoeffPacking);
329  ASSERTL0(localVertexId >= 0 && localVertexId <= 3,
330  "Local Vertex ID must be between 0 and 3");
331  return localVertexId;
332 }
333 
335 {
336  unsigned int i;
337  const unsigned int nBndryCoeff = NumBndryCoeffs();
338 
339  if (outarray.size() != nBndryCoeff)
340  {
341  outarray = Array<OneD, unsigned int>(nBndryCoeff);
342  }
343 
344  for (i = 0; i < nBndryCoeff; i++)
345  {
346  outarray[i] = i;
347  }
348 }
349 
351 {
352  unsigned int i;
353  const unsigned int nBndryCoeff = NumBndryCoeffs();
354 
355  if (outarray.size() != m_ncoeffs - nBndryCoeff)
356  {
357  outarray = Array<OneD, unsigned int>(m_ncoeffs - nBndryCoeff);
358  }
359 
360  for (i = nBndryCoeff; i < m_ncoeffs; i++)
361  {
362  outarray[i - nBndryCoeff] = i;
363  }
364 }
365 
366 //---------------------------------------
367 // Wrapper functions
368 //---------------------------------------
369 
371 {
372  DNekMatSharedPtr Mat;
373 
374  switch (mkey.GetMatrixType())
375  {
376  case eNBasisTrans:
377  Mat = GenNBasisTransMatrix();
378  break;
379  default:
381  break;
382  }
383 
384  return Mat;
385 }
386 
388 {
389  return StdNodalTetExp::v_GenMatrix(mkey);
390 }
391 } // namespace StdRegions
392 } // 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
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void NodalToModalTranspose(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_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool mult=true)
void NodalToModal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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 DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
void GetNodalPoints(Array< OneD, const NekDouble > &x, Array< OneD, const NekDouble > &y, Array< OneD, const NekDouble > &z)
virtual int v_GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
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)
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
LibUtilities::PointsKey m_nodalPointsKey
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:309
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTetExp.cpp:521
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:866
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:675
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:64
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